SU_25_River_Rockfish_Proj_Ambient_V0.3_Fec-Class_Summary

Setup, and Packages

Code Packages:

# Load packages
library(gridExtra) # For using grid.table to save tables as images

## Warning: package 'gridExtra' was built under R version 4.4.3

library(flextable) # Alternative method to save df as image

## Warning: package 'flextable' was built under R version 4.4.3

library(tidyverse)

## Warning: package 'ggplot2' was built under R version 4.4.3

## Warning: package 'lubridate' was built under R version 4.4.3

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.2     ✔ tibble    3.2.1
## ✔ lubridate 1.9.4     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::combine() masks gridExtra::combine()
## ✖ purrr::compose() masks flextable::compose()
## ✖ dplyr::filter()  masks stats::filter()
## ✖ dplyr::lag()     masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(zoo)

## 
## Attaching package: 'zoo'
## 
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

library(lme4)

## Loading required package: Matrix
## 
## Attaching package: 'Matrix'
## 
## The following objects are masked from 'package:tidyr':
## 
##     expand, pack, unpack

library(lmerTest)

## 
## Attaching package: 'lmerTest'
## 
## The following object is masked from 'package:lme4':
## 
##     lmer
## 
## The following object is masked from 'package:stats':
## 
##     step

library(rlang)

## 
## Attaching package: 'rlang'
## 
## The following objects are masked from 'package:purrr':
## 
##     %@%, flatten, flatten_chr, flatten_dbl, flatten_int, flatten_lgl,
##     flatten_raw, invoke, splice

library(nortest) # For ad.test
library(car) # For homogeneity of variance with leveneTest

## Warning: package 'car' was built under R version 4.4.3

## Loading required package: carData

## Warning: package 'carData' was built under R version 4.4.3

## 
## Attaching package: 'car'
## 
## The following object is masked from 'package:dplyr':
## 
##     recode
## 
## The following object is masked from 'package:purrr':
## 
##     some

library(rcompanion) # For scheirerRayHare test

## Warning: package 'rcompanion' was built under R version 4.4.3

library(pastecs) # For stat.desc

## Warning: package 'pastecs' was built under R version 4.4.3

## 
## Attaching package: 'pastecs'
## 
## The following objects are masked from 'package:dplyr':
## 
##     first, last
## 
## The following object is masked from 'package:tidyr':
## 
##     extract

library(ggpubr) # For assumption plots

## Warning: package 'ggpubr' was built under R version 4.4.3

## 
## Attaching package: 'ggpubr'
## 
## The following objects are masked from 'package:flextable':
## 
##     border, font, rotate

library(conover.test) # For the conover.Iman() test function
library(DescTools) # Contain alternative ConoverTest function

## Warning: package 'DescTools' was built under R version 4.4.3

## 
## Attaching package: 'DescTools'
## 
## The following object is masked from 'package:car':
## 
##     Recode

library(rstatix) # tidyverse adapted stat tests and mshapiro_test() functions

## Warning: package 'rstatix' was built under R version 4.4.3

## 
## Attaching package: 'rstatix'
## 
## The following object is masked from 'package:stats':
## 
##     filter

library(DHARMa) # For multi-level lm models and shapiro test

## Warning: package 'DHARMa' was built under R version 4.4.3

## This is DHARMa 0.4.7. For overview type '?DHARMa'. For recent changes, type news(package = 'DHARMa')

library(patchwork) # For combining ggplots

## Warning: package 'patchwork' was built under R version 4.4.3

Data Entry: Read-in data

# Read and call data into df

# Primary samples df
ABL90 <- read.csv("Raw_Data/ABL90_Raw.csv")

# River's Parturition metadata
partuition_subcat <- read.csv("Raw_Data/River's_Rockfish_Metadata_Parturition_V7.csv")

Data Filtering

# Data sifting: ABL90 dataset

# Step 1: Remove missing info
ABL_set1 <- ABL90 %>%
  filter(Patient.ID_edited != "") 
  

# Step 2: Separate Patient.ID into columns of sample types: blood plasma and instant freeze plasma 
ABL_set2 <- separate(ABL_set1, 'Patient.ID_edited', into = c("Patient.ID_edited", "Sample_Type"), sep = ",")  

## Warning: Expected 2 pieces. Missing pieces filled with `NA` in 4 rows [489, 619, 1034,
## 1177].

# Step 3:
# Convert Patient.ID_edited column data to character type
ABL_set2$Patient.ID_edited <- as.character(ABL_set2$Patient.ID_edited) 

# Step 4 & 5: Sussing out specific sample errors  

# Remove insufficient samples
ABL_set3 <- ABL_set2 %>%
  filter(!is.na("Type")) %>%
  filter(!str_detect(Errors.detected.during.measurement, "Insufficient sample"))

# Remove inhomogeneous samples
ABL_set3 <- ABL_set3 %>%
  filter(!is.na("Type")) %>%
  filter(!str_detect(Errors.detected.during.measurement, "Inhomogeneous sample"))

# Step 6: Filter for only blood samples
ABL_b_samp <- ABL_set3 %>%
  filter(Sample_Type == "b")

Merge Datasets:

Join ABL90 df with parturition_subcat df: Making ABL_merged

# Rename columns to match: Change 'ID' in parturition_subcat to 'Patient.ID_edited' so it is ready to join with ABL90 df

# Please note the 'Treatment' col in parturition_subcat does not match with the finalized 'Ambient.Or.OAH' col in metadata_atresia_guide or ABL90_merged df

# Rename 'ID' col to 'patient.ID_edited'
partuition_subcat <- partuition_subcat %>%
  rename(Patient.ID_edited = ID)

# Connect main ABL90 dataset with my (River's) sub categorized parturition metadata 
# ABL_merged <- partuition_subcat %>%
# left_join(ABL_b_samp, by = "Patient.ID_edited")

ABL_merged <- ABL_b_samp %>%
inner_join(partuition_subcat, by = "Patient.ID_edited")

Rename Columns Neatly:

# Rename Treatment & Parturition Outcome
ABL_merged <- ABL_merged %>%
  rename(Ambient.Or.OAH = Consensus_General_Treatment,
         Pregnant.Or.Atresia = Consensus_Atresia_Or_Pregnant)

Remove Samples: Mortalities, NA’s, or Missing Info, and Replicates

Remove Replicate ID’s

# Cross Check all Columns: Cross validation to suss correct replicate row
check_params <- ABL_merged %>%
  select(Patient.ID_edited, Time, Sample.., Measuring.Mode, Ambient.Or.OAH, Pregnant.Or.Atresia, pH, Glu..mmol.L., Hct...., Na...mmol.L., Cl...mmol.L., K...mmol.L., Ca....mmol.L., pCO2..mmHg.) %>%
  filter(Patient.ID_edited == "9782D")
print(check_params)

##   Patient.ID_edited             Time Sample.. Measuring.Mode Ambient.Or.OAH
## 1             9782D 12/28/2023 10:04     1016         C 65uL        Ambient
## 2             9782D   1/9/2024 17:21      863         S 65uL        Ambient
##   Pregnant.Or.Atresia    pH Glu..mmol.L. Hct.... Na...mmol.L. Cl...mmol.L.
## 1    Post Parturition 7.348          0.5     9.8          184          149
## 2    Post Parturition 7.446          0.7    13.3          181          161
##   K...mmol.L. Ca....mmol.L. pCO2..mmHg.
## 1         2.8          1.32         4.2
## 2         2.9          1.31         6.8

check_params <- ABL_merged %>%
  select(Patient.ID_edited, Time, Sample.., Measuring.Mode, Ambient.Or.OAH, Pregnant.Or.Atresia, pH, Glu..mmol.L., Hct...., Na...mmol.L., Cl...mmol.L., K...mmol.L., Ca....mmol.L., pCO2..mmHg.) %>%
  filter(Patient.ID_edited == "9783D")
print(check_params)

##   Patient.ID_edited            Time Sample.. Measuring.Mode Ambient.Or.OAH
## 1             9783D 1/19/2024 10:37      934         S 65uL        Ambient
## 2             9783D  1/7/2024 12:32      860         S 65uL        Ambient
##   Pregnant.Or.Atresia    pH Glu..mmol.L. Hct.... Na...mmol.L. Cl...mmol.L.
## 1    Post Parturition 7.484          1.7     0.2          175          193
## 2    Post Parturition 7.450          1.3    14.7          179          164
##   K...mmol.L. Ca....mmol.L. pCO2..mmHg.
## 1        34.5          1.29         8.5
## 2         2.5          1.30         5.3

# Removed replicates: 9782D (2x) and 9783D (2x)
ABL_merged <- ABL_merged %>%
  filter(Sample.. != "863",
         Sample.. != "934")

Review Replicates:

Row	ID	Time	Sample..	Measuring.Mode	pH	Glu.mmol.L.	Hct….	Na…mmol.L.	Cl…mmol.L.	K…mmol.L.	Ca.mmol.L.
32	9782D	12/28/2023 10:04	1016	C 65uL	7.348	0.5	9.8	184	149	2.8	1.32
39	9782D	1/9/2024 17:21	863	S 65uL	7.446	0.7	13.3	181	161	2.9	1.31
36	9783D	1/19/2024 10:37	934	S 65uL	7.484	1.7	0.2	175	193	34.5	1.29
41	9783D	1/7/2024 12:32	860	S 65uL	7.450	1.3	14.7	179	164	2.5	1.30

Methods for replicate sample removal: Check processing date of IDs, and remove samples that do not match that date.

Case ID_9782D:

• Origianl ABL90 set contains:

• 9782D,b (Sample.. 1016, Time = 12/28/2023 10:04:00 AM)

• 9782D,p (Sample.. 1020, Time = 12/28/2023 10:25:00 AM)

• 9782D,p (Sample.. 865, Time = 1/9/2024 5:29:00 PM)

• 9782D,b (Sample.. 863, Time = 1/9/2024 5:21:00 PM)

• 2024_Fish_Metadata.csv says ID_9782D was processed 3/29/2024

• Note Processing date in metadata will not match time date stamp in ABL90 data due to the machine being reset, and innacurrate dates being automatically recorded.

Deciding info: - For ID-9782D, refer to metadata sheet and look for other samples with that same processing date, then view thoes samples ABL90 ‘Time’. Try to find a matching ‘Time’ with other IDs and 9782D that were processed together during the same time recorded by the machine.

• ID_97838,b (single obs, Sample.. 1024) Time = 12/28/2023 11:49:00 AM

• Since 9782D,b Sample.. 1016 and 97838,b Sample.. 1024 both have a date timestamp occurrence on 12/28/2023, I have decided to keep 9782D,b Sample.. 1016 and drop 9783D,b Sample.. 863

Case ID_9783D

• Original ABL90 dataset contains:

• 9783D,b (Sample.. 934, Time = 1/19/2024 10:37:00 AM)

• 9783D,c (Sample.. 861, Time = 1/7/2024 12:38:00 PM)

• 9783D,b (Sample.. 860, Time = 1/7/2024 12:32:00 PM)

• 2024_Fish_Metadata.csv says ID_9783D was processed 3/20/2024

• Informed that ID_9783D: Sample.. 934 is a plasma sample (on account of its low hct value), remove this and keep Sample.. 860 (parameter values appear correct).

Samples: Data subsets

Remove Motalities and No info IDs and assign analysis ready data-frames:

# New df with Moralities removed: Note none of these samples made it into ABL90 df anyway, so looks like they are already filtered out.
 Live_Samples <- ABL_merged %>%
  filter(Patient.ID_edited != "9780C", # Mortality
         Patient.ID_edited != "777AE", # Mortality
         Patient.ID_edited != "777CA") # Mortality after parturition


# New df with mortality and 'No info' Id's removed
 Primary_Samples <- Live_Samples %>%
   filter(Patient.ID_edited != "777A0", # No info
           Patient.ID_edited != "9782F", # No info
           Patient.ID_edited != "777B3", # No info
           Patient.ID_edited != "777AA", # No info
           Patient.ID_edited != "777DE", # No info
           Patient.ID_edited != "777CE", # No info
           Patient.ID_edited != "777A6") # No info
 
 
# New df of Only Ambient Treatment: For testing parturition success
 Ambient_Samples <- Primary_Samples %>%
   filter(Ambient.Or.OAH == "Ambient")
 
# New df of Fecundity Samples:
 Fecundity_Samples <- Primary_Samples %>%
   filter(Fecundity_Or_Brood_Count != "NA",
     Fecundity_Class != "NA") # removes 97706

# Fecundity samples without atresia Ids
 Fecundity_No_Atresia_Samples <- Primary_Samples %>%
   filter(All_Fecundity != "Na",
          All_Fecundity != 0)  # removes atresia IDs

# Fecundity Ambient Samples
 Amb_Fec_Samples <- Primary_Samples %>%
   filter(Ambient.Or.OAH == "Ambient",
          All_Fecundity != "NA") # removes OAH treatment samples and ID 97706
 
# Fecundity Ambient Samples without atresia Ids
 Amb_Fec_No_Atresia_Samples <- Primary_Samples %>%
   filter(Ambient.Or.OAH == "Ambient", # removes OAH treatment 
          All_Fecundity != "NA", # removes all missing info IDs
          All_Fecundity != 0) # removes atresia IDs

Change Data Types:

# Change Continuous Data to type to numeric, double, or factor

# Change data type of ions to factor or numeric
#glimpse(General_Samples)
#glimpse(Ambient_Samples)
 
Primary_Samples <- Primary_Samples %>%
  mutate(across(c(All_Fecundity, pH, Hct...., Glu..mmol.L., Na...mmol.L., Cl...mmol.L., K...mmol.L., Ca....mmol.L.), as.numeric))

Ambient_Samples <- Ambient_Samples %>%
  mutate(across(c(All_Fecundity, pH, Hct...., Glu..mmol.L., Na...mmol.L., Cl...mmol.L., K...mmol.L., Ca....mmol.L.), as.numeric))

Fecundity_Samples <- Fecundity_Samples %>%
  mutate(across(c(All_Fecundity, pH, Hct...., Glu..mmol.L., Na...mmol.L., Cl...mmol.L., K...mmol.L., Ca....mmol.L.), as.numeric))

Fecundity_No_Atresia_Samples <- Fecundity_Samples %>%
  mutate(across(c(All_Fecundity, pH, Hct...., Glu..mmol.L., Na...mmol.L., Cl...mmol.L., K...mmol.L., Ca....mmol.L.), as.numeric))

Amb_Fec_Samples <- Amb_Fec_Samples %>%
  mutate(across(c(All_Fecundity, pH, Hct...., Glu..mmol.L., Na...mmol.L., Cl...mmol.L., K...mmol.L., Ca....mmol.L.), as.numeric))
  
Amb_Fec_No_Atresia_Samples <- Amb_Fec_No_Atresia_Samples %>%
  mutate(across(c(All_Fecundity, pH, Hct...., Glu..mmol.L., Na...mmol.L., Cl...mmol.L., K...mmol.L., Ca....mmol.L.), as.numeric))

# glimpse(Primary_Samples)
# glimpse(Ambient_Samples)
# glimpse(Fecundity_Samples)
# glimpse(Amb_Fec_Samples)

Order Factor Levels

# Fecundity data 

#unique(Fecundity_Samples$Fecundity_Class)

# Arrange the order of parturition categories for visualizations & tests

# Primary Samples

# Arrange Treatment
Primary_Samples$Ambient.Or.OAH <- ordered(Primary_Samples$Ambient.Or.OAH, levels = c("Ambient", "OAH"))

# Arrange Parturition outcome
Primary_Samples$Pregnant.Or.Atresia <- ordered(Primary_Samples$Pregnant.Or.Atresia, levels = c("Post Parturition", "Atresia"))

# Arrange Brood Condition
Primary_Samples$Consensus_Brood_Condition <- ordered(Primary_Samples$Consensus_Brood_Condition, levels = c("Excellent", "Normal", "Low", "Very Low", "Atresia"))

# Arrange Fecundity Class
Primary_Samples$Fecundity_Class <- ordered(Primary_Samples$Fecundity_Class, levels = c("High (>50,000)", "Low (~1,000s)", "Atresia"))

# Remove NAs from ordered character factors
Primary_Samples <- Primary_Samples %>%
  filter(Fecundity_Class != is.na(Fecundity_Class)) %>%
  filter(Consensus_Brood_Condition != is.na(Consensus_Brood_Condition))
 
 
# Ambient data
Ambient_Samples$Ambient.Or.OAH <- ordered(Ambient_Samples$Ambient.Or.OAH, levels = c("Ambient", "OAH"))

Ambient_Samples$Pregnant.Or.Atresia <- ordered(Ambient_Samples$Pregnant.Or.Atresia, levels = c("Post Parturition", "Atresia"))

Ambient_Samples$Consensus_Brood_Condition <- ordered(Ambient_Samples$Consensus_Brood_Condition, levels = c("Excellent", "Normal", "Low", "Atresia"))

Ambient_Samples$Fecundity_Class <- ordered(Ambient_Samples$Fecundity_Class, levels = c("High (>50,000)", "Low (~1,000s)", "Atresia"))

# Remove NAs from ordered character factors
Ambient_Samples <- Ambient_Samples %>%
  filter(Fecundity_Class != is.na(Fecundity_Class)) %>%
  filter(Consensus_Brood_Condition != is.na(Consensus_Brood_Condition))

# Ambient Fecundity data
Amb_Fec_No_Atresia_Samples$Ambient.Or.OAH <- ordered(Amb_Fec_No_Atresia_Samples$Ambient.Or.OAH, levels = c("Ambient", "OAH"))

Amb_Fec_No_Atresia_Samples$Pregnant.Or.Atresia <- ordered(Amb_Fec_No_Atresia_Samples$Pregnant.Or.Atresia, levels = c("Post Parturition", "Atresia"))

Amb_Fec_No_Atresia_Samples$Consensus_Brood_Condition <- ordered(Amb_Fec_No_Atresia_Samples$Consensus_Brood_Condition, levels = c("Excellent", "Normal", "Low", "Atresia"))

Amb_Fec_No_Atresia_Samples$Fecundity_Class <- ordered(Amb_Fec_No_Atresia_Samples$Fecundity_Class, levels = c("High (>50,000)", "Low (~1,000s)", "Atresia"))

# Remove NAs from ordered character factors
Amb_Fec_No_Atresia_Samples <- Amb_Fec_No_Atresia_Samples %>%
  filter(Fecundity_Class != is.na(Fecundity_Class)) %>%
  filter(Consensus_Brood_Condition != is.na(Consensus_Brood_Condition))

Save filtered dataset into data-worked folder

# Save merged dataset in a worked folder

write.csv(x = Ambient_Samples, file = "Worked-Data/Ambient_Samples")

Sample Size: n

# Sample Size

# Ambient data:
ambient_sample_size_table <- Ambient_Samples %>%
  group_by(Ambient.Or.OAH, Consensus_Brood_Condition) %>%
  summarise(n_size = n_distinct(Patient.ID_edited), .groups = "drop")
print(ambient_sample_size_table)

## # A tibble: 4 × 3
##   Ambient.Or.OAH Consensus_Brood_Condition n_size
##   <ord>          <ord>                      <int>
## 1 Ambient        Excellent                      1
## 2 Ambient        Normal                         4
## 3 Ambient        Low                            3
## 4 Ambient        Atresia                       13

ambient_fecundity_class_sample_size_table <- Ambient_Samples %>%
  group_by(Ambient.Or.OAH, Fecundity_Class) %>%
  summarise(n_size = n_distinct(Patient.ID_edited), .groups = "drop")
print(ambient_fecundity_class_sample_size_table)

## # A tibble: 3 × 3
##   Ambient.Or.OAH Fecundity_Class n_size
##   <ord>          <ord>            <int>
## 1 Ambient        High (>50,000)       4
## 2 Ambient        Low (~1,000s)        4
## 3 Ambient        Atresia             13

amb_fec_no_atresia_sample_size_table <- Amb_Fec_No_Atresia_Samples %>%
  group_by(Ambient.Or.OAH, Weight_Adjusted_Fecundity_Fecundity.Mean_Weight) %>%
  summarise(n_size = n_distinct(Patient.ID_edited), .groups = "drop")
print(amb_fec_no_atresia_sample_size_table)

## # A tibble: 8 × 3
##   Ambient.Or.OAH Weight_Adjusted_Fecundity_Fecundity.Mean_Weight n_size
##   <ord>                                                    <dbl>  <int>
## 1 Ambient                                                   6.22      1
## 2 Ambient                                                  68.8       1
## 3 Ambient                                                 114.        1
## 4 Ambient                                                 114.        1
## 5 Ambient                                                 148.        1
## 6 Ambient                                                 161.        1
## 7 Ambient                                                 204.        1
## 8 Ambient                                                 356.        1

Sample Size Barplot:

# View Sample Distributions
# Sample Size barplot
# Ambient Samples

# Consensus brood condition sample barplot
ambient_sample_size_barplot <- Ambient_Samples %>%
  group_by(Consensus_Brood_Condition, Ambient.Or.OAH) %>% 
  summarise(n_size = n_distinct(Patient.ID_edited), .groups = "drop") %>%
  ggplot(aes(x = Consensus_Brood_Condition, y = n_size)) +
  geom_bar(stat = "identity") +
  facet_grid(~ Ambient.Or.OAH) +
  scale_y_continuous(expand = c(0, 0), limits = c(0, 15.5), breaks = c(0, 3, 6, 9, 12, 15)) +
  labs(title = "Ambient Sample Size",
        x = "Parturition Outcome",
        y = "Sample Size (n)") +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        title = element_text(size = 12),
        plot.title = element_text(size = 30, hjust = 0.5),
        axis.text.x = element_text(size = 10),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))
print(ambient_sample_size_barplot)

# Fecundity class sample barplot
ambient_sample_size_barplot2 <- Ambient_Samples %>%
  group_by(Fecundity_Class, Ambient.Or.OAH) %>% 
  summarise(n_size = n_distinct(Patient.ID_edited), .groups = "drop") %>%
  ggplot(aes(x = Fecundity_Class, y = n_size)) +
  geom_bar(stat = "identity") +
  facet_grid(~ Ambient.Or.OAH) +
  scale_y_continuous(expand = c(0, 0), limits = c(0, 15.5), breaks = c(0, 3, 6, 9, 12, 15)) +
  labs(title = "Ambient Sample Size",
        x = "Parturition Outcome",
        y = "Sample Size (n)") +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        title = element_text(size = 12),
        plot.title = element_text(size = 30, hjust = 0.5),
        axis.text.x = element_text(size = 10),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))
print(ambient_sample_size_barplot2)

# Amb Fec No atresia sample size

amb_fec_no_atresia_barplot <- Amb_Fec_No_Atresia_Samples %>%
  group_by(Fecundity_Class, Ambient.Or.OAH) %>% 
  summarise(n_size = n_distinct(Patient.ID_edited), .groups = "drop") %>%
  ggplot(aes(x = Fecundity_Class, y = n_size)) +
  geom_bar(stat = "identity") +
  facet_grid(~ Ambient.Or.OAH) +
  scale_y_continuous(expand = c(0, 0), limits = c(0, 15.5), breaks = c(4, 8, 12)) +
  labs(title = "Ambient No Atresia Sample Size",
        x = "Parturition Outcome",
        y = "Sample Size (n)") +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        title = element_text(size = 12),
        plot.title = element_text(size = 30, hjust = 0.5),
        axis.text.x = element_text(size = 10),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))
print(amb_fec_no_atresia_barplot)

Paramater Analysis:

Co2 : Carbon Dioxide

CO2 Boxplot

pCO2_ambient_fec_class_boxplot <- ggplot(data = Ambient_Samples) + 
  geom_boxplot(aes(x = Fecundity_Class, y = pCO2..mmHg.), fill = "#189bff") +
  geom_point(aes(x = Fecundity_Class, y = pCO2..mmHg.), alpha = 0.5, colour = "black") +
  labs(title = "Carbon Dioxide", x = "Fecundity Class", y = "pCO2 (mmHg)") +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        title = element_text(size = 12),
        plot.title = element_text(size = 30, hjust = 0.5),
        axis.text.x = element_text(size = 10),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_blank(),
        axis.title.y = element_text(size = 20))

print(pCO2_ambient_fec_class_boxplot)

# Ambient No Atresia Samples Scatterplot
pCO2_amb_fec_no_atresia_scatterplot <- Amb_Fec_No_Atresia_Samples %>%
  ggplot() + 
  geom_point(aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = pCO2..mmHg.), shape = 21, fill = "#189bff") +
  labs(title = "Carbon Dioxide", x = "Weight Adjusted Fecundity", y = "pCO2 (mmHg)") +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        title = element_text(size = 12),
        plot.title = element_text(size = 30, hjust = 0.5),
        axis.text.x = element_text(size = 20),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))

print(pCO2_amb_fec_no_atresia_scatterplot)

pCo2 Analysis

Differences: Between Fecundity Class

# pCO2
# Ambient ANOVA Test

# Testing for differences in ambient parturition success metrics
# One-way ANOVA (parametric) or Kruskal-Wallis test (nonparametric)

pCO2_aov_ambient_fec <- aov(pCO2..mmHg. ~ Fecundity_Class, data = Ambient_Samples)
summary(pCO2_aov_ambient_fec)

##                 Df Sum Sq Mean Sq F value Pr(>F)
## Fecundity_Class  2  10.10   5.050   1.817  0.191
## Residuals       18  50.02   2.779

# If significant, Post Hoc Test: 
# - Tukey's HSD (parametric, normality & equal variance across groups)
# - Dunn's test (non-parametric, after Kruskal-Wallis test)

TukeyHSD(pCO2_aov_ambient_fec)

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = pCO2..mmHg. ~ Fecundity_Class, data = Ambient_Samples)
## 
## $Fecundity_Class
##                                    diff        lwr      upr     p adj
## Low (~1,000s)-High (>50,000) -1.2750000 -4.2833026 1.733303 0.5371486
## Atresia-High (>50,000)        0.5384615 -1.8940746 2.970998 0.8401814
## Atresia-Low (~1,000s)         1.8134615 -0.6190746 4.245998 0.1666659

# Install flexplot
# devtools::install_github("dustinfife/flexplot", ref="development")
require(flexplot)

## Loading required package: flexplot

## 
## Attaching package: 'flexplot'

## The following object is masked from 'package:ggplot2':
## 
##     flip_data

pCO2_ambient_fec_sim <- simulateResiduals(fittedModel = pCO2_aov_ambient_fec) 

## Warning in checkModel(fittedModel): DHARMa: fittedModel not in class of
## supported models. Absolutely no guarantee that this will work!

plot(pCO2_ambient_fec_sim)

plotResiduals(pCO2_ambient_fec_sim)

testDispersion(pCO2_ambient_fec_sim)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.92221, p-value = 0.888
## alternative hypothesis: two.sided

testOutliers(pCO2_ambient_fec_sim)  

## 
##  DHARMa outlier test based on exact binomial test with approximate
##  expectations
## 
## data:  pCO2_ambient_fec_sim
## outliers at both margin(s) = 0, observations = 21, p-value = 1
## alternative hypothesis: true probability of success is not equal to 0.007968127
## 95 percent confidence interval:
##  0.0000000 0.1610976
## sample estimates:
## frequency of outliers (expected: 0.00796812749003984 ) 
##                                                      0

outlierTest(pCO2_aov_ambient_fec)

## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
##   rstudent unadjusted p-value Bonferroni p
## 7 2.080913           0.052875           NA

# Parametric tests

# Normality
shapiro.test(Ambient_Samples$pCO2..mmHg.)

## 
##  Shapiro-Wilk normality test
## 
## data:  Ambient_Samples$pCO2..mmHg.
## W = 0.97964, p-value = 0.9201

shapiro.test(residuals(aov(pCO2..mmHg. ~ Fecundity_Class, data = Ambient_Samples)))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(aov(pCO2..mmHg. ~ Fecundity_Class, data = Ambient_Samples))
## W = 0.95328, p-value = 0.3922

# QQplot:
pCO2_ambient_fec_res_qqplot <- ggqqplot(residuals(aov(pCO2..mmHg. ~ Fecundity_Class, data = Ambient_Samples))) +
labs(title = "Ambient pCO2 Residual QQplot",
     subtitle = "residuals(aov(pCO2 ~ Fecundity_Class, data = Ambient_Samples))",
                               x = "pCO2 Theoretical Quantiles (Predicted values)",
                               y = "pCO2 Sample Quantiles") +
  theme(plot.title = element_text(hjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5))
print(pCO2_ambient_fec_res_qqplot)

# Parametric Variance test (Bartlett's test): require at least 2 obs for each group
bartlett.test(pCO2..mmHg. ~ Fecundity_Class, data = Ambient_Samples)

## 
##  Bartlett test of homogeneity of variances
## 
## data:  pCO2..mmHg. by Fecundity_Class
## Bartlett's K-squared = 0.88695, df = 2, p-value = 0.6418

# Nonparametric variance test (Levene's test):
leveneTest(pCO2..mmHg. ~ Fecundity_Class, data = Ambient_Samples)

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  2  0.5105 0.6086
##       18

# Nonparametric stat test: Kruskall Wallis
kruskal.test(pCO2..mmHg. ~ Fecundity_Class, data = Ambient_Samples)

## 
##  Kruskal-Wallis rank sum test
## 
## data:  pCO2..mmHg. by Fecundity_Class
## Kruskal-Wallis chi-squared = 2.6912, df = 2, p-value = 0.2604

# Nonparametric Post Hoc: Dunn's test
DunnTest(pCO2..mmHg. ~ Fecundity_Class, data = Ambient_Samples)

## 
##  Dunn's test of multiple comparisons using rank sums : holm  
## 
##                              mean.rank.diff   pval    
## Low (~1,000s)-High (>50,000)      -4.125000 0.6921    
## Atresia-High (>50,000)             1.673077 0.6921    
## Atresia-Low (~1,000s)              5.798077 0.3043    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlations: Parameter to Numerical Fecundity (without Atresia IDs)

# pCO2
# Linear Model
pCO2_amb_fec_no_atresia_lm <- lm(pCO2..mmHg. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = Ambient_Samples)
summary(pCO2_amb_fec_no_atresia_lm)

## 
## Call:
## lm(formula = pCO2..mmHg. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, 
##     data = Ambient_Samples)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -2.760 -1.069  0.120  1.401  1.880 
## 
## Coefficients:
##                                                 Estimate Std. Error t value
## (Intercept)                                     4.770067   1.138478   4.190
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight 0.001314   0.006485   0.203
##                                                 Pr(>|t|)   
## (Intercept)                                      0.00575 **
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight  0.84608   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.777 on 6 degrees of freedom
##   (13 observations deleted due to missingness)
## Multiple R-squared:  0.006801,   Adjusted R-squared:  -0.1587 
## F-statistic: 0.04108 on 1 and 6 DF,  p-value: 0.8461

pCO2_ambient_sim <- simulateResiduals(fittedModel = pCO2_amb_fec_no_atresia_lm) 
plot(pCO2_ambient_sim)

plotResiduals(pCO2_ambient_sim)

testDispersion(pCO2_ambient_sim)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.84221, p-value = 0.92
## alternative hypothesis: two.sided

testOutliers(pCO2_ambient_sim)  

## 
##  DHARMa outlier test based on exact binomial test with approximate
##  expectations
## 
## data:  pCO2_ambient_sim
## outliers at both margin(s) = 0, observations = 8, p-value = 1
## alternative hypothesis: true probability of success is not equal to 0.007968127
## 95 percent confidence interval:
##  0.0000000 0.3694166
## sample estimates:
## frequency of outliers (expected: 0.00796812749003984 ) 
##                                                      0

outlierTest(pCO2_amb_fec_no_atresia_lm)

## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
##     rstudent unadjusted p-value Bonferroni p
## 21 -2.262971           0.073079      0.58463

# Normality test
shapiro.test(residuals(pCO2_amb_fec_no_atresia_lm))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(pCO2_amb_fec_no_atresia_lm)
## W = 0.93831, p-value = 0.5945

# qq plot
ggqqplot(residuals(lm(pCO2..mmHg. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = Ambient_Samples)))

Residuals Plot (if significant)

# pCO2
# Ambient samples

# ggplot with geom_smooth(aes(x=,y=)), method = "lm")

# Original lm model (no transformation)
pCO2_ambient_lm_scatterplot <- ggplot(data = Ambient_Samples, aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = pCO2..mmHg.)) +
  geom_point() +
  geom_smooth(aes(group = 1), method = "lm") +
  labs(title = "Carbon Dioxide", x = "Weight Adjusted Fecundity", y = "pCO2 (mmHg)", color = "Category") +
  scale_y_continuous() +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        plot.title = element_text(size = 25, hjust = 0.5),
        axis.text.x = element_text(size = 15),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))

print(pCO2_ambient_lm_scatterplot)

## `geom_smooth()` using formula = 'y ~ x'

## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).

## Warning: Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).

pH

pH boxplot

pH_ambient_fec_class_boxplot <- ggplot(data = Ambient_Samples) + 
  geom_boxplot(aes(x = Fecundity_Class, y = pH), fill = "#189bff") +
  geom_point(aes(x = Fecundity_Class, y = pH), alpha = 0.5, colour = "black") +
  labs(title = "pH", x = "Fecundity Class", y = "pH") +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        title = element_text(size = 12),
        plot.title = element_text(size = 30, hjust = 0.5),
        axis.text.x = element_text(size = 10),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_blank(),
        axis.title.y = element_text(size = 20))

print(pH_ambient_fec_class_boxplot)

# Ambient No Atresia Samples Scatterplot
pH_amb_fec_no_atresia_scatterplot <- Amb_Fec_No_Atresia_Samples %>%
  ggplot() + 
  geom_point(aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = pH), shape = 21, fill = "#189bff") +
  labs(title = "pH", x = "Weight Adjusted Fecundity", y = "pH") +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        title = element_text(size = 12),
        plot.title = element_text(size = 30, hjust = 0.5),
        axis.text.x = element_text(size = 20),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))

print(pH_amb_fec_no_atresia_scatterplot)

pH Stat Tests

Differences: Between Fecundity Class

# pH
# Ambient ANOVA Test

# Testing for differences in ambient parturition success metrics
# One-way ANOVA (parametric) or Kruskal-Wallis test (nonparametric)

pH_aov_ambient_fec <- aov(pH ~ Fecundity_Class, data = Ambient_Samples)
summary(pH_aov_ambient_fec)

##                 Df  Sum Sq  Mean Sq F value Pr(>F)
## Fecundity_Class  2 0.02183 0.010916   1.336  0.288
## Residuals       18 0.14711 0.008173

# If significant, Post Hoc Test: 
# - Tukey's HSD (parametric, normality & equal variance across groups)
# - Dunn's test (non-parametric, after Kruskal-Wallis test)

TukeyHSD(pH_aov_ambient_fec)

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = pH ~ Fecundity_Class, data = Ambient_Samples)
## 
## $Fecundity_Class
##                                     diff        lwr        upr     p adj
## Low (~1,000s)-High (>50,000)  0.02675000 -0.1363959 0.18989587 0.9084726
## Atresia-High (>50,000)       -0.05080769 -0.1827287 0.08111329 0.5966404
## Atresia-Low (~1,000s)        -0.07755769 -0.2094787 0.05436329 0.3141611

pH_ambient_fec_sim <- simulateResiduals(fittedModel = pH_aov_ambient_fec) 

## Warning in checkModel(fittedModel): DHARMa: fittedModel not in class of
## supported models. Absolutely no guarantee that this will work!

plot(pH_ambient_fec_sim)

plotResiduals(pH_ambient_fec_sim)

testDispersion(pH_ambient_fec_sim)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.91709, p-value = 0.872
## alternative hypothesis: two.sided

testOutliers(pH_ambient_fec_sim)  

## 
##  DHARMa outlier test based on exact binomial test with approximate
##  expectations
## 
## data:  pH_ambient_fec_sim
## outliers at both margin(s) = 0, observations = 21, p-value = 1
## alternative hypothesis: true probability of success is not equal to 0.007968127
## 95 percent confidence interval:
##  0.0000000 0.1610976
## sample estimates:
## frequency of outliers (expected: 0.00796812749003984 ) 
##                                                      0

outlierTest(pH_aov_ambient_fec)

## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
##    rstudent unadjusted p-value Bonferroni p
## 13  3.11884          0.0062471      0.13119

# Parametric tests

# Normality
shapiro.test(Ambient_Samples$pH)

## 
##  Shapiro-Wilk normality test
## 
## data:  Ambient_Samples$pH
## W = 0.98534, p-value = 0.9808

shapiro.test(residuals(aov(pH ~ Fecundity_Class, data = Ambient_Samples)))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(aov(pH ~ Fecundity_Class, data = Ambient_Samples))
## W = 0.97088, p-value = 0.7524

# QQplot:
pH_ambient_fec_res_qqplot <- ggqqplot(residuals(aov(pH ~ Fecundity_Class, data = Ambient_Samples))) +
labs(title = "Ambient pH Residual QQplot",
     subtitle = "residuals(aov(pH ~ Fecundity_Class, data = Ambient_Samples))",
                               x = "pH Theoretical Quantiles (Predicted values)",
                               y = "pH Sample Quantiles") +
  theme(plot.title = element_text(hjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5))
print(pH_ambient_fec_res_qqplot)

# Parametric Variance test (Bartlett's test): require at least 2 obs for each group
bartlett.test(pH ~ Fecundity_Class, data = Ambient_Samples)

## 
##  Bartlett test of homogeneity of variances
## 
## data:  pH by Fecundity_Class
## Bartlett's K-squared = 1.4521, df = 2, p-value = 0.4838

# Nonparametric variance test (Levene's test):
leveneTest(pH ~ Fecundity_Class, data = Ambient_Samples)

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  2    0.44 0.6508
##       18

# Nonparametric stat test: Kruskall Wallis
kruskal.test(pH ~ Fecundity_Class, data = Ambient_Samples)

## 
##  Kruskal-Wallis rank sum test
## 
## data:  pH by Fecundity_Class
## Kruskal-Wallis chi-squared = 3.0375, df = 2, p-value = 0.219

# Nonparametric Post Hoc: Dunn's test
DunnTest(pH ~ Fecundity_Class, data = Ambient_Samples)

## 
##  Dunn's test of multiple comparisons using rank sums : holm  
## 
##                              mean.rank.diff   pval    
## Low (~1,000s)-High (>50,000)       2.250000 0.6424    
## Atresia-High (>50,000)            -3.519231 0.6424    
## Atresia-Low (~1,000s)             -5.769231 0.3117    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlations: Parameter to Numerical Fecundity (without Atresia IDs)

# Ph
# Linear Model
pH_amb_fec_no_atresia_lm <- lm(pH ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = Ambient_Samples)
summary(pH_amb_fec_no_atresia_lm)

## 
## Call:
## lm(formula = pH ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, 
##     data = Ambient_Samples)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.11361 -0.02694  0.01591  0.04102  0.05887 
## 
## Coefficients:
##                                                   Estimate Std. Error t value
## (Intercept)                                      7.4734874  0.0403735 185.109
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight -0.0002228  0.0002300  -0.969
##                                                 Pr(>|t|)    
## (Intercept)                                     1.68e-12 ***
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight     0.37    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.06302 on 6 degrees of freedom
##   (13 observations deleted due to missingness)
## Multiple R-squared:  0.1352, Adjusted R-squared:  -0.008899 
## F-statistic: 0.9383 on 1 and 6 DF,  p-value: 0.3701

pH_ambient_sim <- simulateResiduals(fittedModel = pH_amb_fec_no_atresia_lm) 
plot(pH_ambient_sim)

## qu = 0.25, log(sigma) = -2.642568 : outer Newton did not converge fully.

plotResiduals(pH_amb_fec_no_atresia_lm)

## qu = 0.25, log(sigma) = -2.642568 : outer Newton did not converge fully.

testDispersion(pH_amb_fec_no_atresia_lm)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.84221, p-value = 0.92
## alternative hypothesis: two.sided

testOutliers(pH_ambient_sim)  

## 
##  DHARMa outlier test based on exact binomial test with approximate
##  expectations
## 
## data:  pH_ambient_sim
## outliers at both margin(s) = 0, observations = 8, p-value = 1
## alternative hypothesis: true probability of success is not equal to 0.007968127
## 95 percent confidence interval:
##  0.0000000 0.3694166
## sample estimates:
## frequency of outliers (expected: 0.00796812749003984 ) 
##                                                      0

outlierTest(pH_amb_fec_no_atresia_lm)

## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
##    rstudent unadjusted p-value Bonferroni p
## 3 -2.862653           0.035297      0.28237

# Normality test
shapiro.test(residuals(pH_amb_fec_no_atresia_lm))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(pH_amb_fec_no_atresia_lm)
## W = 0.90496, p-value = 0.32

# qq plot
ggqqplot(residuals(lm(pH ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = Ambient_Samples)))

# Transform data: square root (undo log scale?)
pH_ambient_sqrt <- Amb_Fec_No_Atresia_Samples %>%
  mutate(pH = sqrt(pH))

# Rerun lm model with square root transformed data
pH_amb_fec_no_atresia_lm_sqrt <- lm(pH ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = pH_ambient_sqrt)
summary(pH_amb_fec_no_atresia_lm_sqrt)

## 
## Call:
## lm(formula = pH ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, 
##     data = pH_ambient_sqrt)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.020889 -0.004920  0.002934  0.007526  0.010787 
## 
## Coefficients:
##                                                   Estimate Std. Error t value
## (Intercept)                                      2.734e+00  7.412e-03 368.821
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight -4.078e-05  4.222e-05  -0.966
##                                                 Pr(>|t|)    
## (Intercept)                                     2.68e-14 ***
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight    0.371    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.01157 on 6 degrees of freedom
## Multiple R-squared:  0.1346, Adjusted R-squared:  -0.009642 
## F-statistic: 0.9331 on 1 and 6 DF,  p-value: 0.3714

pH_amb_fec_sqrt_sim <- simulateResiduals(fittedModel = pH_amb_fec_no_atresia_lm_sqrt) 
plot(pH_amb_fec_sqrt_sim)

## qu = 0.25, log(sigma) = -2.642951 : outer Newton did not converge fully.

plotResiduals(pH_amb_fec_no_atresia_lm_sqrt)

## qu = 0.25, log(sigma) = -2.642951 : outer Newton did not converge fully.

testDispersion(pH_amb_fec_no_atresia_lm_sqrt)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.84221, p-value = 0.92
## alternative hypothesis: two.sided

# Residuals of model for shapiro test
pH_amb_fec_no_atresia_lm_res <- residuals(pH_amb_fec_no_atresia_lm_sqrt) 

# Normality test
shapiro.test(pH_amb_fec_no_atresia_lm_res) # normal

## 
##  Shapiro-Wilk normality test
## 
## data:  pH_amb_fec_no_atresia_lm_res
## W = 0.90414, p-value = 0.3147

# Parametric variance test: bartlett's test
# Require at least 2 obs for each group
#bartlett.test(pH ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = Fecundity_No_Atresia_Samples) 

# Nonparametric variance test (Levene's test)
# NA for quantitative explanatory variables
# leveneTest(pH ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = Fecundity_No_Atresia_Samples)

Residuals Plot (if significant)

# pH
# Ambient samples

# ggplot with geom_smooth(aes(x=,y=)), method = "lm")

# Original lm model (no transformation)
pH_ambient_lm_scatterplot <- ggplot(data = Ambient_Samples, aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = pH)) +
  geom_point() +
  geom_smooth(aes(group = 1), method = "lm") +
  labs(title = "pH", x = "Weight Adjusted Fecundity", y = "pH", color = "Category") +
  scale_y_continuous() +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        plot.title = element_text(size = 25, hjust = 0.5),
        axis.text.x = element_text(size = 15),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))

print(pH_ambient_lm_scatterplot)

## `geom_smooth()` using formula = 'y ~ x'

## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).

## Warning: Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).

Hematocrit

# hct

# Hct boxplot
hct_ambient_fec_class_boxplot <- ggplot(data = Ambient_Samples) + 
  geom_boxplot(aes(x = Fecundity_Class, y = Hct....), fill = "#189bff") +
  geom_point(aes(x = Fecundity_Class, y = Hct....), alpha = 0.5, colour = "black") +
  labs(title = "Hematocrit", x = "Fecundity", y = "Hematocrit (%)") +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        title = element_text(size = 12),
        plot.title = element_text(size = 30, hjust = 0.5),
        axis.text.x = element_text(size = 10),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_blank(),
        axis.title.y = element_text(size = 20))

print(hct_ambient_fec_class_boxplot)

# Ambient No Atresia Sample Scatterplot
hct_amb_fec_no_atresia_scatterplot <- Amb_Fec_No_Atresia_Samples %>%
  ggplot() + 
  geom_point(aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = Hct....), shape = 21, fill = "#189bff") +
  labs(title = "Hematocrit", x = "Weight Adjusted Fecundity", y = "Hematocrit (%)") +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        title = element_text(size = 12),
        plot.title = element_text(size = 30, hjust = 0.5),
        axis.text.x = element_text(size = 20),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))

print(hct_amb_fec_no_atresia_scatterplot)

Hct Stat Tests

Differences: Between Fecundity Class

# hct
# Ambient ANOVA Test

# Testing for differences in ambient parturition success metrics
# One-way ANOVA (parametric) or Kruskal-Wallis test (nonparametric)

hct_aov_ambient_fec <- aov(Hct.... ~ Fecundity_Class, data = Ambient_Samples)
summary(hct_aov_ambient_fec)

##                 Df Sum Sq Mean Sq F value Pr(>F)  
## Fecundity_Class  2  73.83   36.92   4.824  0.021 *
## Residuals       18 137.74    7.65                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# If significant, Post Hoc Test: 
# - Tukey's HSD (parametric, normality & equal variance across groups)
# - Dunn's test (non-parametric, after Kruskal-Wallis test)

TukeyHSD(hct_aov_ambient_fec)

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Hct.... ~ Fecundity_Class, data = Ambient_Samples)
## 
## $Fecundity_Class
##                              diff        lwr      upr     p adj
## Low (~1,000s)-High (>50,000) -3.4 -8.3921454 1.592145 0.2186627
## Atresia-High (>50,000)        1.5 -2.5366864 5.536686 0.6176987
## Atresia-Low (~1,000s)         4.9  0.8633136 8.936686 0.0162654

hct_ambient_fec_sim <- simulateResiduals(fittedModel = hct_aov_ambient_fec) 

## Warning in checkModel(fittedModel): DHARMa: fittedModel not in class of
## supported models. Absolutely no guarantee that this will work!

plot(hct_ambient_fec_sim)

plotResiduals(hct_ambient_fec_sim)

testDispersion(hct_ambient_fec_sim)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.91709, p-value = 0.872
## alternative hypothesis: two.sided

testOutliers(hct_ambient_fec_sim)  

## 
##  DHARMa outlier test based on exact binomial test with approximate
##  expectations
## 
## data:  hct_ambient_fec_sim
## outliers at both margin(s) = 0, observations = 21, p-value = 1
## alternative hypothesis: true probability of success is not equal to 0.007968127
## 95 percent confidence interval:
##  0.0000000 0.1610976
## sample estimates:
## frequency of outliers (expected: 0.00796812749003984 ) 
##                                                      0

outlierTest(hct_aov_ambient_fec)

## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
##    rstudent unadjusted p-value Bonferroni p
## 20 2.329698           0.032404      0.68048

# Parametric tests

# Normality
shapiro.test(Ambient_Samples$Hct....)

## 
##  Shapiro-Wilk normality test
## 
## data:  Ambient_Samples$Hct....
## W = 0.97482, p-value = 0.8355

shapiro.test(residuals(aov(Hct.... ~ Fecundity_Class, data = Ambient_Samples)))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(aov(Hct.... ~ Fecundity_Class, data = Ambient_Samples))
## W = 0.97094, p-value = 0.7538

# QQplot:
hct_ambient_fec_res_qqplot <- ggqqplot(residuals(aov(Hct.... ~ Fecundity_Class, data = Ambient_Samples))) +
labs(title = "Ambient Hematocrit Residual QQplot",
     subtitle = "residuals(aov(Hct.... ~ Fecundity_Class, data = Ambient_Samples))",
                               x = "Hct Theoretical Quantiles (Predicted values)",
                               y = "Hct Sample Quantiles") +
  theme(plot.title = element_text(hjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5))
print(hct_ambient_fec_res_qqplot)

# Parametric Variance test (Bartlett's test): require at least 2 obs for each group
bartlett.test(Hct.... ~ Fecundity_Class, data = Ambient_Samples)

## 
##  Bartlett test of homogeneity of variances
## 
## data:  Hct.... by Fecundity_Class
## Bartlett's K-squared = 1.7219, df = 2, p-value = 0.4228

# Nonparametric variance test (Levene's test):
leveneTest(Hct.... ~ Fecundity_Class, data = Ambient_Samples)

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  2  1.2003 0.3241
##       18

# Nonparametric stat test: Kruskall Wallis
kruskal.test(Hct.... ~ Fecundity_Class, data = Ambient_Samples)

## 
##  Kruskal-Wallis rank sum test
## 
## data:  Hct.... by Fecundity_Class
## Kruskal-Wallis chi-squared = 6.7821, df = 2, p-value = 0.03367

# Nonparametric Post Hoc: Dunn's Test
DunnTest(Hct.... ~ Fecundity_Class, data = Ambient_Samples)

## 
##  Dunn's test of multiple comparisons using rank sums : holm  
## 
##                              mean.rank.diff   pval    
## Low (~1,000s)-High (>50,000)      -6.750000 0.2476    
## Atresia-High (>50,000)             2.480769 0.4843    
## Atresia-Low (~1,000s)              9.230769 0.0277 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlations: Parameter to Numerical Fecundity (without Atresia IDs)

# hct
# Linear Model
hct_amb_fec_no_atresia_lm <- lm(Hct.... ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = Ambient_Samples)
summary(hct_amb_fec_no_atresia_lm)

## 
## Call:
## lm(formula = Hct.... ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, 
##     data = Ambient_Samples)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.5587 -1.5142 -0.7795  0.5244  7.1866 
## 
## Coefficients:
##                                                 Estimate Std. Error t value
## (Intercept)                                     12.81211    2.32735   5.505
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight  0.01494    0.01326   1.127
##                                                 Pr(>|t|)   
## (Intercept)                                      0.00151 **
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight  0.30267   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.633 on 6 degrees of freedom
##   (13 observations deleted due to missingness)
## Multiple R-squared:  0.1748, Adjusted R-squared:  0.03725 
## F-statistic: 1.271 on 1 and 6 DF,  p-value: 0.3027

hct_ambient_sim <- simulateResiduals(fittedModel = hct_amb_fec_no_atresia_lm) 
plot(hct_ambient_sim)

## 
##  DHARMa: qgam was unable to calculate quantile regression for quantile 0.25. Possibly to few (unique) data points / predictions. The quantile will be ommited in plots and significance calculations.

## We had to increase `err` for some of the quantiles. See fit$calibr$err

plotResiduals(hct_amb_fec_no_atresia_lm)

## 
##  DHARMa: qgam was unable to calculate quantile regression for quantile 0.25. Possibly to few (unique) data points / predictions. The quantile will be ommited in plots and significance calculations. 
## 
## We had to increase `err` for some of the quantiles. See fit$calibr$err

testDispersion(hct_amb_fec_no_atresia_lm)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.84221, p-value = 0.92
## alternative hypothesis: two.sided

testOutliers(hct_ambient_sim)  

## 
##  DHARMa outlier test based on exact binomial test with approximate
##  expectations
## 
## data:  hct_ambient_sim
## outliers at both margin(s) = 0, observations = 8, p-value = 1
## alternative hypothesis: true probability of success is not equal to 0.007968127
## 95 percent confidence interval:
##  0.0000000 0.3694166
## sample estimates:
## frequency of outliers (expected: 0.00796812749003984 ) 
##                                                      0

outlierTest(hct_amb_fec_no_atresia_lm)

## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
##    rstudent unadjusted p-value Bonferroni p
## 20 3.953321           0.010815     0.086518

# Normality test
shapiro.test(residuals(hct_amb_fec_no_atresia_lm))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(hct_amb_fec_no_atresia_lm)
## W = 0.86076, p-value = 0.1222

# qq plot
ggqqplot(residuals(lm(Hct.... ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = Ambient_Samples)))

# Transform data: square root (undo log scale?)
hct_ambient_sqrt <- Ambient_Samples %>%
  mutate(Hct.... = sqrt(Hct....))

# Rerun lm model with square root transformed data
hct_amb_fec_no_atresia_lm_sqrt <- lm(Hct.... ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = hct_ambient_sqrt)
summary(hct_amb_fec_no_atresia_lm_sqrt)

## 
## Call:
## lm(formula = Hct.... ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, 
##     data = hct_ambient_sqrt)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.45751 -0.18537 -0.07738  0.06130  0.87682 
## 
## Coefficients:
##                                                 Estimate Std. Error t value
## (Intercept)                                     3.549914   0.292111  12.153
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight 0.002034   0.001664   1.223
##                                                 Pr(>|t|)    
## (Intercept)                                     1.89e-05 ***
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight    0.267    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.456 on 6 degrees of freedom
##   (13 observations deleted due to missingness)
## Multiple R-squared:  0.1995, Adjusted R-squared:  0.06603 
## F-statistic: 1.495 on 1 and 6 DF,  p-value: 0.2673

hct_amb_fec_sqrt_sim <- simulateResiduals(fittedModel = hct_amb_fec_no_atresia_lm_sqrt) 
plot(hct_amb_fec_sqrt_sim)

plotResiduals(hct_amb_fec_no_atresia_lm_sqrt)

testDispersion(hct_amb_fec_no_atresia_lm_sqrt)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.84221, p-value = 0.92
## alternative hypothesis: two.sided

# Residuals of model for shapiro test
hct_amb_sqrt_res <- residuals(hct_amb_fec_no_atresia_lm_sqrt) 

# Normality test
shapiro.test(hct_amb_sqrt_res) 

## 
##  Shapiro-Wilk normality test
## 
## data:  hct_amb_sqrt_res
## W = 0.87306, p-value = 0.1615

# Parametric variance test: bartlett's test
# Require at least 2 obs for each group
#bartlett.test(Hct.... ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = hct_ambient_sqrt) 

# Nonparametric variance test (Levene's test)
# NA for quantitative explanatory variables
# leveneTest(Hct.... ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = hct_ambient_sqrt)

Residuals Plot (if significant)

# hct
# Ambient samples

# ggplot with geom_smooth(aes(x=,y=)), method = "lm")

hct_ambient_lm_scatterplot <- ggplot(data = Ambient_Samples, aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = Hct....)) +
  geom_point() +
  geom_smooth(aes(group = 1), method = "lm") +
  labs(title = "Hematocrit", x = "Weight Adjusted Fecundity", y = "Hematocrit (%)", color = "Category") +
  scale_y_continuous() +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        plot.title = element_text(size = 25, hjust = 0.5),
        axis.text.x = element_text(size = 15),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))

print(hct_ambient_lm_scatterplot)

## `geom_smooth()` using formula = 'y ~ x'

## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).

## Warning: Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).

Glucose

Glucose graphs

# Glucose

# Glucose boxplot 
glucose_ambient_fec_class_boxplot <- ggplot(data = Ambient_Samples) + 
  geom_boxplot(aes(x = Fecundity_Class, y = Glu..mmol.L.), fill = "#189bff") +
  geom_point(aes(x = Fecundity_Class, y = Glu..mmol.L.), alpha = 0.5, colour = "black") +
  labs(title = "Glucose", x = "Fecundity Class", y = "Glucose (mmol/L)") +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        title = element_text(size = 12),
        plot.title = element_text(size = 30, hjust = 0.5),
        axis.text.x = element_text(size = 10),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_blank(),
        axis.title.y = element_text(size = 20))

print(glucose_ambient_fec_class_boxplot)

# Glucose Ambient No Atresia Scatterplot
glucose_amb_fec_no_atresia_scatterplot <- Amb_Fec_No_Atresia_Samples %>%
  ggplot() + 
  geom_point(aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = Glu..mmol.L.), shape = 21, fill = "#189bff") +
  labs(title = "Glucose", x = "Weight Adjusted Fecundity", y = "Glucose (mmol/L)") +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        title = element_text(size = 12),
        plot.title = element_text(size = 30, hjust = 0.5),
        axis.text.x = element_text(size = 20),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))

print(glucose_amb_fec_no_atresia_scatterplot)

Glucose Stat Tests

Differences: Between Fecundity Class

# glucose
# Ambient ANOVA Test

# Testing for differences in ambient parturition success metrics
# One-way ANOVA (parametric) or Kruskal-Wallis test (nonparametric)

glucose_aov_ambient_fec <- aov(Glu..mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
summary(glucose_aov_ambient_fec)

##                 Df Sum Sq Mean Sq F value Pr(>F)
## Fecundity_Class  2  1.208  0.6041   1.317  0.293
## Residuals       18  8.257  0.4587

# If significant, Post Hoc Test: 
# - Tukey's HSD (parametric, normality & equal variance across groups)
# - Dunn's test (non-parametric, after Kruskal-Wallis test)

TukeyHSD(glucose_aov_ambient_fec)

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Glu..mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
## 
## $Fecundity_Class
##                               diff        lwr      upr     p adj
## Low (~1,000s)-High (>50,000) 0.125 -1.0973103 1.347310 0.9632210
## Atresia-High (>50,000)       0.550 -0.4383693 1.538369 0.3518867
## Atresia-Low (~1,000s)        0.425 -0.5633693 1.413369 0.5277964

glu_ambient_fec_sim <- simulateResiduals(fittedModel = glucose_aov_ambient_fec) 

## Warning in checkModel(fittedModel): DHARMa: fittedModel not in class of
## supported models. Absolutely no guarantee that this will work!

plot(glu_ambient_fec_sim)

plotResiduals(glu_ambient_fec_sim)

testDispersion(glu_ambient_fec_sim)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.91709, p-value = 0.872
## alternative hypothesis: two.sided

testOutliers(glu_ambient_fec_sim)  

## 
##  DHARMa outlier test based on exact binomial test with approximate
##  expectations
## 
## data:  glu_ambient_fec_sim
## outliers at both margin(s) = 0, observations = 21, p-value = 1
## alternative hypothesis: true probability of success is not equal to 0.007968127
## 95 percent confidence interval:
##  0.0000000 0.1610976
## sample estimates:
## frequency of outliers (expected: 0.00796812749003984 ) 
##                                                      0

outlierTest(glucose_aov_ambient_fec)

##    rstudent unadjusted p-value Bonferroni p
## 13 4.830979          0.0001563    0.0032822

# Parametric tests

# Normality
shapiro.test(Ambient_Samples$Glu..mmol.L.)

## 
##  Shapiro-Wilk normality test
## 
## data:  Ambient_Samples$Glu..mmol.L.
## W = 0.85634, p-value = 0.005473

shapiro.test(residuals(aov(Glu..mmol.L. ~ Fecundity_Class, data = Ambient_Samples)))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(aov(Glu..mmol.L. ~ Fecundity_Class, data = Ambient_Samples))
## W = 0.82056, p-value = 0.001376

# QQplot:
glucose_ambient_fec_res_qqplot <- ggqqplot(residuals(aov(Glu..mmol.L. ~ Fecundity_Class, data = Ambient_Samples))) +
labs(title = "Ambient Glucose Interactive Residual QQplot",
     subtitle = "residuals(aov(Glu..mmol.L. ~ Fecundity_Class, data = Ambient_Samples))",
                               x = "Glucose Theoretical Quantiles (Predicted values)",
                               y = "Glucose Sample Quantiles") +
  theme(plot.title = element_text(hjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5))
print(glucose_ambient_fec_res_qqplot)

# Parametric Variance test (Bartlett's test): require at least 2 obs for each group
bartlett.test(Glu..mmol.L. ~ Fecundity_Class, data = Ambient_Samples)

## 
##  Bartlett test of homogeneity of variances
## 
## data:  Glu..mmol.L. by Fecundity_Class
## Bartlett's K-squared = 3.669, df = 2, p-value = 0.1597

# Nonparametric variance test (Levene's test):
leveneTest(Glu..mmol.L. ~ Fecundity_Class, data = Ambient_Samples)

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  2  0.9391 0.4093
##       18

# Nonparametric stat test: Kruskall Wallis
kruskal.test(Glu..mmol.L. ~ Fecundity_Class, data = Ambient_Samples)

## 
##  Kruskal-Wallis rank sum test
## 
## data:  Glu..mmol.L. by Fecundity_Class
## Kruskal-Wallis chi-squared = 4.0755, df = 2, p-value = 0.1303

# Nonparametric Post Hoc: Dunn's test
DunnTest(Glu..mmol.L. ~ Fecundity_Class, data = Ambient_Samples)

## 
##  Dunn's test of multiple comparisons using rank sums : holm  
## 
##                              mean.rank.diff   pval    
## Low (~1,000s)-High (>50,000)       2.625000 0.5480    
## Atresia-High (>50,000)             6.663462 0.1778    
## Atresia-Low (~1,000s)              4.038462 0.5060    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlations: Parameter to Numerical Fecundity (without Atresia IDs)

# glucose
# Linear Model
glucose_amb_fec_no_atresia_lm <- lm(Glu..mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = Ambient_Samples)
summary(glucose_amb_fec_no_atresia_lm)

## 
## Call:
## lm(formula = Glu..mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, 
##     data = Ambient_Samples)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.71290 -0.36526 -0.05942  0.16524  1.19856 
## 
## Coefficients:
##                                                  Estimate Std. Error t value
## (Intercept)                                     1.5622030  0.4085226   3.824
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight 0.0003435  0.0023269   0.148
##                                                 Pr(>|t|)   
## (Intercept)                                      0.00872 **
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight  0.88746   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6377 on 6 degrees of freedom
##   (13 observations deleted due to missingness)
## Multiple R-squared:  0.00362,    Adjusted R-squared:  -0.1624 
## F-statistic: 0.0218 on 1 and 6 DF,  p-value: 0.8875

glucose_ambient_sim <- simulateResiduals(fittedModel = glucose_amb_fec_no_atresia_lm) 
plot(glucose_ambient_sim)

## Warning in newton(lsp = lsp, X = G$X, y = G$y, Eb = G$Eb, UrS = G$UrS, L = G$L,
## : Fitting terminated with step failure - check results carefully

plotResiduals(glucose_amb_fec_no_atresia_lm)

## Warning in newton(lsp = lsp, X = G$X, y = G$y, Eb = G$Eb, UrS = G$UrS, L = G$L,
## : Fitting terminated with step failure - check results carefully

testDispersion(glucose_amb_fec_no_atresia_lm)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.84221, p-value = 0.92
## alternative hypothesis: two.sided

testOutliers(glucose_ambient_sim)  

## 
##  DHARMa outlier test based on exact binomial test with approximate
##  expectations
## 
## data:  glucose_ambient_sim
## outliers at both margin(s) = 0, observations = 8, p-value = 1
## alternative hypothesis: true probability of success is not equal to 0.007968127
## 95 percent confidence interval:
##  0.0000000 0.3694166
## sample estimates:
## frequency of outliers (expected: 0.00796812749003984 ) 
##                                                      0

outlierTest(glucose_amb_fec_no_atresia_lm)

## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
##   rstudent unadjusted p-value Bonferroni p
## 5 3.287331           0.021778      0.17422

# Normality test:
shapiro.test(residuals(glucose_amb_fec_no_atresia_lm))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(glucose_amb_fec_no_atresia_lm)
## W = 0.92103, p-value = 0.4383

# qq plot
ggqqplot(residuals(lm(Glu..mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = Ambient_Samples)))

# Transform data: square root 
glucose_ambient_sqrt <- Amb_Fec_No_Atresia_Samples %>%
  mutate(Glu..mmol.L. = sqrt(Glu..mmol.L.))

# Rerun lm model with square root transformed data
glucose_amb_fec_no_atresia_lm_sqrt <- lm(Glu..mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = glucose_ambient_sqrt)
summary(glucose_amb_fec_no_atresia_lm_sqrt)

## 
## Call:
## lm(formula = Glu..mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, 
##     data = glucose_ambient_sqrt)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.30380 -0.13647 -0.01130  0.07862  0.42723 
## 
## Coefficients:
##                                                  Estimate Std. Error t value
## (Intercept)                                     1.2241905  0.1551624   7.890
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight 0.0001917  0.0008838   0.217
##                                                 Pr(>|t|)    
## (Intercept)                                      0.00022 ***
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight  0.83544    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2422 on 6 degrees of freedom
## Multiple R-squared:  0.007783,   Adjusted R-squared:  -0.1576 
## F-statistic: 0.04706 on 1 and 6 DF,  p-value: 0.8354

glucose_amb_fec_sqrt_sim <- simulateResiduals(fittedModel = glucose_amb_fec_no_atresia_lm_sqrt) 
plot(glucose_amb_fec_sqrt_sim)

plotResiduals(glucose_amb_fec_no_atresia_lm_sqrt)

testDispersion(glucose_amb_fec_no_atresia_lm_sqrt)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.84221, p-value = 0.92
## alternative hypothesis: two.sided

# Residuals of model for shapiro test
glucose_amb_fec_no_atresia_lm_res <- residuals(glucose_amb_fec_no_atresia_lm_sqrt) 

# Normality test
shapiro.test(glucose_amb_fec_no_atresia_lm_res) 

## 
##  Shapiro-Wilk normality test
## 
## data:  glucose_amb_fec_no_atresia_lm_res
## W = 0.96004, p-value = 0.8104

# Parametric variance test: bartlett's test
# Require at least 2 obs for each group
#bartlett.test(Glu..mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = glucose_ambient_sqrt) 

# Nonparametric variance test (Levene's test)
# NA for quantitative explanatory variables
# leveneTest(Glu..mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = glucose_ambient_sqrt)

Residuals Plot (if significant)

# glucose
# Ambient samples

# ggplot with geom_smooth(aes(x=,y=)), method = "lm")

glucose_ambient_lm_scatterplot <- ggplot(data = Ambient_Samples, aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = Glu..mmol.L.)) +
  geom_point() +
  geom_smooth(aes(group = 1), method = "lm") +
  labs(title = "Glucose", x = "Weight Adjusted Fecundity", y = "Glucose (mmol/L)", color = "Category") +
  scale_y_continuous() +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        plot.title = element_text(size = 25, hjust = 0.5),
        axis.text.x = element_text(size = 15),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))

print(glucose_ambient_lm_scatterplot)

## `geom_smooth()` using formula = 'y ~ x'

## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).

## Warning: Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).

Sodium

# Na+

# sodium boxplot

sodium_ambient_fec_class_boxplot <- ggplot(data = Ambient_Samples) + 
  geom_boxplot(aes(x = Fecundity_Class, y = Na...mmol.L.), fill = "#189bff") +
  geom_point(aes(x = Fecundity_Class, y = Na...mmol.L.), alpha = 0.5, colour = "black") +
  labs(title = "Sodium", x = "Fecundity Class", y = "Sodium (mmol/L)") +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        title = element_text(size = 12),
        plot.title = element_text(size = 30, hjust = 0.5),
        axis.text.x = element_text(size = 10),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_blank(),
        axis.title.y = element_text(size = 20))

print(sodium_ambient_fec_class_boxplot)

# sodium ambient no atresia scatterplot

sodium_amb_fec_no_atresia_scatterplot <- Ambient_Samples %>%
  ggplot() + 
  geom_point(aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = Na...mmol.L.), shape = 21, fill = "#189bff") +
  labs(title = "Sodium", x = "Weight Adjusted Fecundity", y = "Sodium (mmol/L)") +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        plot.title = element_text(size = 30, hjust = 0.5),
        axis.text.x = element_text(size = 20),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))

print(sodium_amb_fec_no_atresia_scatterplot)

## Warning: Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).

Na+ Stat Tests

Differences: Between Fecundity Class

# Na+
# Ambient ANOVA Test

# Testing for differences in ambient parturition success metrics
# One-way ANOVA (parametric) or Kruskal-Wallis test (nonparametric)

sodium_aov_ambient_fec <- aov(Na...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
summary(sodium_aov_ambient_fec)

##                 Df Sum Sq Mean Sq F value Pr(>F)
## Fecundity_Class  2  156.8   78.38   0.865  0.438
## Residuals       18 1631.8   90.66

# If significant, Post Hoc Test: 
# - Tukey's HSD (parametric, normality & equal variance across groups)
# - Dunn's test (non-parametric, after Kruskal-Wallis test)

TukeyHSD(sodium_aov_ambient_fec)

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Na...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
## 
## $Fecundity_Class
##                                   diff       lwr      upr     p adj
## Low (~1,000s)-High (>50,000) -2.000000 -19.18271 15.18271 0.9526465
## Atresia-High (>50,000)        4.480769  -9.41330 18.37484 0.6939179
## Atresia-Low (~1,000s)         6.480769  -7.41330 20.37484 0.4737484

na_ambient_fec_sim <- simulateResiduals(fittedModel = sodium_aov_ambient_fec) 

## Warning in checkModel(fittedModel): DHARMa: fittedModel not in class of
## supported models. Absolutely no guarantee that this will work!

plot(na_ambient_fec_sim)

plotResiduals(na_ambient_fec_sim)

testDispersion(na_ambient_fec_sim)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.91709, p-value = 0.872
## alternative hypothesis: two.sided

testOutliers(na_ambient_fec_sim)  

## 
##  DHARMa outlier test based on exact binomial test with approximate
##  expectations
## 
## data:  na_ambient_fec_sim
## outliers at both margin(s) = 1, observations = 21, p-value = 0.1546
## alternative hypothesis: true probability of success is not equal to 0.007968127
## 95 percent confidence interval:
##  0.001204883 0.238159910
## sample estimates:
## frequency of outliers (expected: 0.00796812749003984 ) 
##                                             0.04761905

outlierTest(sodium_aov_ambient_fec)

##    rstudent unadjusted p-value Bonferroni p
## 17 12.20437         7.7681e-10   1.6313e-08

# Parametric tests

# Normality
shapiro.test(Ambient_Samples$Na...mmol.L.)

## 
##  Shapiro-Wilk normality test
## 
## data:  Ambient_Samples$Na...mmol.L.
## W = 0.52676, p-value = 3.625e-07

shapiro.test(residuals(aov(Na...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(aov(Na...mmol.L. ~ Fecundity_Class, data = Ambient_Samples))
## W = 0.56321, p-value = 8.235e-07

# QQplot:
sodium_ambient_fec_res_qqplot <- ggqqplot(residuals(aov(Na...mmol.L. ~ Fecundity_Class, data = Ambient_Samples))) +
labs(title = "Ambient Sodium Residual QQplot",
     subtitle = "residuals(aov(Na...mmol.L. ~ Fecundity_Class, data = Ambient_Samples))",
                               x = "Sodium Theoretical Quantiles (Predicted values)",
                               y = "Sodium Sample Quantiles") +
  theme(plot.title = element_text(hjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5))
print(sodium_ambient_fec_res_qqplot)

# Parametric Variance test (Bartlett's test): require at least 2 obs for each group
bartlett.test(Na...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)

## 
##  Bartlett test of homogeneity of variances
## 
## data:  Na...mmol.L. by Fecundity_Class
## Bartlett's K-squared = 7.4634, df = 2, p-value = 0.02395

# Nonparametric variance test (Levene's test):
leveneTest(Na...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  2  0.2098 0.8127
##       18

# Nonparametric stat test: Kruskall Wallis
kruskal.test(Na...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)

## 
##  Kruskal-Wallis rank sum test
## 
## data:  Na...mmol.L. by Fecundity_Class
## Kruskal-Wallis chi-squared = 2.8608, df = 2, p-value = 0.2392

# Nonparametric Post Hoc: Dunn's test
DunnTest(Na...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)

## 
##  Dunn's test of multiple comparisons using rank sums : holm  
## 
##                              mean.rank.diff   pval    
## Low (~1,000s)-High (>50,000)      -3.375000 0.8730    
## Atresia-High (>50,000)             2.451923 0.8730    
## Atresia-Low (~1,000s)              5.826923 0.2899    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlations: Parameter to Numerical Fecundity (without Atresia IDs)

# Na+
# Linear Model
sodium_amb_fec_no_atresia_lm <- lm(Na...mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = Ambient_Samples)
summary(sodium_amb_fec_no_atresia_lm)

## 
## Call:
## lm(formula = Na...mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, 
##     data = Ambient_Samples)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.7526 -1.0177  0.5473  2.5230  3.4276 
## 
## Coefficients:
##                                                  Estimate Std. Error t value
## (Intercept)                                     1.754e+02  2.355e+00  74.474
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight 2.289e-03  1.342e-02   0.171
##                                                 Pr(>|t|)    
## (Intercept)                                     3.94e-10 ***
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight     0.87    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.677 on 6 degrees of freedom
##   (13 observations deleted due to missingness)
## Multiple R-squared:  0.004827,   Adjusted R-squared:  -0.161 
## F-statistic: 0.0291 on 1 and 6 DF,  p-value: 0.8701

sodium_ambient_sim <- simulateResiduals(fittedModel = sodium_amb_fec_no_atresia_lm) 
plot(sodium_ambient_sim)

## 
##  DHARMa: qgam was unable to calculate quantile regression for quantile 0.25. Possibly to few (unique) data points / predictions. The quantile will be ommited in plots and significance calculations.

## 
##  DHARMa: qgam was unable to calculate quantile regression for quantile 0.5. Possibly to few (unique) data points / predictions. The quantile will be ommited in plots and significance calculations.

## 
##  DHARMa: qgam was unable to calculate quantile regression for quantile 0.75. Possibly to few (unique) data points / predictions. The quantile will be ommited in plots and significance calculations.

plotResiduals(sodium_ambient_sim)

## 
##  DHARMa: qgam was unable to calculate quantile regression for quantile 0.25. Possibly to few (unique) data points / predictions. The quantile will be ommited in plots and significance calculations.

## 
##  DHARMa: qgam was unable to calculate quantile regression for quantile 0.5. Possibly to few (unique) data points / predictions. The quantile will be ommited in plots and significance calculations.

## 
##  DHARMa: qgam was unable to calculate quantile regression for quantile 0.75. Possibly to few (unique) data points / predictions. The quantile will be ommited in plots and significance calculations.

testDispersion(sodium_ambient_sim)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.84221, p-value = 0.92
## alternative hypothesis: two.sided

testOutliers(sodium_ambient_sim)  

## 
##  DHARMa outlier test based on exact binomial test with approximate
##  expectations
## 
## data:  sodium_ambient_sim
## outliers at both margin(s) = 0, observations = 8, p-value = 1
## alternative hypothesis: true probability of success is not equal to 0.007968127
## 95 percent confidence interval:
##  0.0000000 0.3694166
## sample estimates:
## frequency of outliers (expected: 0.00796812749003984 ) 
##                                                      0

outlierTest(sodium_amb_fec_no_atresia_lm)

## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
##    rstudent unadjusted p-value Bonferroni p
## 2 -2.997851           0.030174      0.24139

# Normailty test
shapiro.test(residuals(sodium_amb_fec_no_atresia_lm))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(sodium_amb_fec_no_atresia_lm)
## W = 0.89794, p-value = 0.2768

# qq plot
ggqqplot(residuals(lm(Na...mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = Ambient_Samples)))

# Transform data: square root 
sodium_ambient_sqrt <- Amb_Fec_No_Atresia_Samples %>%
  mutate(Na...mmol.L. = sqrt(Na...mmol.L.))

# Rerun lm model with square root transformed data
sodium_amb_fec_no_atresia_lm_sqrt <- lm(Na...mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = sodium_ambient_sqrt)
summary(sodium_amb_fec_no_atresia_lm_sqrt)

## 
## Call:
## lm(formula = Na...mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, 
##     data = sodium_ambient_sqrt)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.25662 -0.03796  0.02115  0.09537  0.12926 
## 
## Coefficients:
##                                                  Estimate Std. Error t value
## (Intercept)                                     1.324e+01  8.923e-02  148.42
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight 8.630e-05  5.083e-04    0.17
##                                                 Pr(>|t|)    
## (Intercept)                                     6.31e-12 ***
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight    0.871    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1393 on 6 degrees of freedom
## Multiple R-squared:  0.004782,   Adjusted R-squared:  -0.1611 
## F-statistic: 0.02883 on 1 and 6 DF,  p-value: 0.8708

sodium_amb_fec_sqrt_sim <- simulateResiduals(fittedModel = sodium_amb_fec_no_atresia_lm_sqrt) 
plot(sodium_amb_fec_sqrt_sim)

## 
##  DHARMa: qgam was unable to calculate quantile regression for quantile 0.25. Possibly to few (unique) data points / predictions. The quantile will be ommited in plots and significance calculations.

## 
##  DHARMa: qgam was unable to calculate quantile regression for quantile 0.5. Possibly to few (unique) data points / predictions. The quantile will be ommited in plots and significance calculations.

## 
##  DHARMa: qgam was unable to calculate quantile regression for quantile 0.75. Possibly to few (unique) data points / predictions. The quantile will be ommited in plots and significance calculations.

plotResiduals(sodium_amb_fec_no_atresia_lm_sqrt)

## 
##  DHARMa: qgam was unable to calculate quantile regression for quantile 0.25. Possibly to few (unique) data points / predictions. The quantile will be ommited in plots and significance calculations.

## 
##  DHARMa: qgam was unable to calculate quantile regression for quantile 0.5. Possibly to few (unique) data points / predictions. The quantile will be ommited in plots and significance calculations.

## 
##  DHARMa: qgam was unable to calculate quantile regression for quantile 0.75. Possibly to few (unique) data points / predictions. The quantile will be ommited in plots and significance calculations.

testDispersion(sodium_amb_fec_no_atresia_lm_sqrt)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.84221, p-value = 0.92
## alternative hypothesis: two.sided

# Residuals of model for shapiro test
sodium_amb_fec_no_atresia_lm_res <- residuals(sodium_amb_fec_no_atresia_lm_sqrt) 

# Normality test
shapiro.test(sodium_amb_fec_no_atresia_lm_res) 

## 
##  Shapiro-Wilk normality test
## 
## data:  sodium_amb_fec_no_atresia_lm_res
## W = 0.89593, p-value = 0.2654

# Parametric variance test: bartlett's test
# Require at least 2 obs for each group
#bartlett.test(Na...mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = sodium_ambient_sqrt) 

# Nonparametric variance test (Levene's test)
# NA for quantitative explanatory variables
# leveneTest(Na...mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = sodium_ambient_sqrt)

Residuals Plot (if significant)

# Na+
# Ambient samples

# ggplot with geom_smooth(aes(x=,y=)), method = "lm")

sodium_ambient_lm_scatterplot <- ggplot(data = Ambient_Samples, aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = Na...mmol.L.)) +
  geom_point() +
  geom_smooth(aes(group = 1), method = "lm") +
  labs(title = "Sodium", x = "Weight Adjusted Fecundity", y = "Sodium (mmol/L)", color = "Category") +
  scale_y_continuous() +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        plot.title = element_text(size = 25, hjust = 0.5),
        axis.text.x = element_text(size = 15),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))

print(sodium_ambient_lm_scatterplot)

## `geom_smooth()` using formula = 'y ~ x'

## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).

## Warning: Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).

Chloride

Chloride graphs

# Cl- 

# Chloride boxplot 
chloride_ambient_fec_class_boxplot <- ggplot(data = Ambient_Samples) + 
  geom_boxplot(aes(x = Fecundity_Class, y = Cl...mmol.L.), fill = "#189bff") +
  geom_point(aes(x = Fecundity_Class, y = Cl...mmol.L.), alpha = 0.5, colour = "black") +
  labs(title = "Chloride", x = "Fecundity Class", y = "Chloride (mmol/L)") +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        title = element_text(size = 12),
        plot.title = element_text(size = 30, hjust = 0.5),
        axis.text.x = element_text(size = 10),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_blank(),
        axis.title.y = element_text(size = 20))

print(chloride_ambient_fec_class_boxplot)

# Chloride Ambient No Atresia Scatterplot
chloride_amb_fec_no_atresia_scatterplot <- Amb_Fec_No_Atresia_Samples %>%
  ggplot() + 
  geom_point(aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = Na...mmol.L.), shape = 21, fill = "#189bff") +
  labs(title = "Chloride", x = "Weight Adjusted Fecundity", y = "Chloride (mmol/L)") +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        plot.title = element_text(size = 30, hjust = 0.5),
        axis.text.x = element_text(size = 20),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))

print(chloride_amb_fec_no_atresia_scatterplot)

Cl- Stat Tests

Differences: Between Fecundity Class

# Cl-
# Ambient ANOVA Test

# Testing for differences in ambient parturition success metrics
# One-way ANOVA (parametric) or Kruskal-Wallis test (nonparametric)

chloride_aov_ambient_fec <- aov(Cl...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
summary(chloride_aov_ambient_fec)

##                 Df Sum Sq Mean Sq F value Pr(>F)
## Fecundity_Class  2    3.0   1.504   0.084   0.92
## Residuals       18  322.2  17.902

# If significant, Post Hoc Test: 
# - Tukey's HSD (parametric, normality & equal variance across groups)
# - Dunn's test (non-parametric, after Kruskal-Wallis test)

TukeyHSD(chloride_aov_ambient_fec)

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Cl...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
## 
## $Fecundity_Class
##                                    diff       lwr      upr     p adj
## Low (~1,000s)-High (>50,000) -0.5000000 -8.135556 7.135556 0.9847331
## Atresia-High (>50,000)        0.4615385 -5.712630 6.635707 0.9801564
## Atresia-Low (~1,000s)         0.9615385 -5.212630 7.135707 0.9170009

cl_ambient_fec_sim <- simulateResiduals(fittedModel = chloride_aov_ambient_fec) 

## Warning in checkModel(fittedModel): DHARMa: fittedModel not in class of
## supported models. Absolutely no guarantee that this will work!

plot(cl_ambient_fec_sim)

plotResiduals(cl_ambient_fec_sim)

testDispersion(cl_ambient_fec_sim)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.91709, p-value = 0.872
## alternative hypothesis: two.sided

testOutliers(cl_ambient_fec_sim)  

## 
##  DHARMa outlier test based on exact binomial test with approximate
##  expectations
## 
## data:  cl_ambient_fec_sim
## outliers at both margin(s) = 0, observations = 21, p-value = 1
## alternative hypothesis: true probability of success is not equal to 0.007968127
## 95 percent confidence interval:
##  0.0000000 0.1610976
## sample estimates:
## frequency of outliers (expected: 0.00796812749003984 ) 
##                                                      0

outlierTest(chloride_aov_ambient_fec)

## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
##    rstudent unadjusted p-value Bonferroni p
## 16 2.003759           0.061299           NA

# Parametric tests

# Normality
shapiro.test(Ambient_Samples$Cl...mmol.L.)

## 
##  Shapiro-Wilk normality test
## 
## data:  Ambient_Samples$Cl...mmol.L.
## W = 0.96502, p-value = 0.6224

shapiro.test(residuals(aov(Cl...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(aov(Cl...mmol.L. ~ Fecundity_Class, data = Ambient_Samples))
## W = 0.96555, p-value = 0.6339

# QQplot:
chloride_ambient_fec_res_qqplot <- ggqqplot(residuals(aov(Cl...mmol.L. ~ Fecundity_Class, data = Ambient_Samples))) +
labs(title = "Ambient Chloride Residual QQplot",
     subtitle = "residuals(aov(Cl...mmol.L. ~ Fecundity_Class, data = Ambient_Samples))",
                               x = "Chloride Theoretical Quantiles (Predicted values)",
                               y = "Chloride Sample Quantiles") +
  theme(plot.title = element_text(hjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5))
print(chloride_ambient_fec_res_qqplot)

# Parametric Variance test (Bartlett's test): require at least 2 obs for each group
bartlett.test(Cl...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)

## 
##  Bartlett test of homogeneity of variances
## 
## data:  Cl...mmol.L. by Fecundity_Class
## Bartlett's K-squared = 0.82247, df = 2, p-value = 0.6628

# Nonparametric variance test (Levene's test):
leveneTest(Cl...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  2  0.5546 0.5838
##       18

# Nonparametric stat test: Kruskall Wallis
kruskal.test(Cl...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)

## 
##  Kruskal-Wallis rank sum test
## 
## data:  Cl...mmol.L. by Fecundity_Class
## Kruskal-Wallis chi-squared = 0.4005, df = 2, p-value = 0.8185

# Nonparametric Post Hoc: Dunns
DunnTest(Cl...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)

## 
##  Dunn's test of multiple comparisons using rank sums : holm  
## 
##                              mean.rank.diff   pval    
## Low (~1,000s)-High (>50,000)     -2.0000000 1.0000    
## Atresia-High (>50,000)            0.2115385 1.0000    
## Atresia-Low (~1,000s)             2.2115385 1.0000    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Remove Outlier

# Cl-

# Test outlier

cl_ambient_fec_sim <- simulateResiduals(fittedModel = chloride_aov_ambient_fec) 

## Warning in checkModel(fittedModel): DHARMa: fittedModel not in class of
## supported models. Absolutely no guarantee that this will work!

plot(cl_ambient_fec_sim)

plotResiduals(cl_ambient_fec_sim)

testDispersion(cl_ambient_fec_sim)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.91709, p-value = 0.872
## alternative hypothesis: two.sided

testOutliers(cl_ambient_fec_sim)  

## 
##  DHARMa outlier test based on exact binomial test with approximate
##  expectations
## 
## data:  cl_ambient_fec_sim
## outliers at both margin(s) = 0, observations = 21, p-value = 1
## alternative hypothesis: true probability of success is not equal to 0.007968127
## 95 percent confidence interval:
##  0.0000000 0.1610976
## sample estimates:
## frequency of outliers (expected: 0.00796812749003984 ) 
##                                                      0

outlierTest(chloride_aov_ambient_fec)

## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
##    rstudent unadjusted p-value Bonferroni p
## 16 2.003759           0.061299           NA

# Find Outlier
boxplot(Ambient_Samples$Cl...mmol.L., horizontal = TRUE)

ggplot(data = Ambient_Samples, aes(x = Fecundity_Class, y = Cl...mmol.L., fill = Ambient.Or.OAH)) +
  geom_boxplot(aes(x = Fecundity_Class, y = Cl...mmol.L.)) +
  geom_point(aes(color = Ambient.Or.OAH), position = position_dodge(width = 0.75), alpha = 0.3, color = "black")

# Remove Outlier 
# cl_amb_no_outlier <- Ambient_Samples %>%
#   filter(Cl...mmol.L. !=)

# Rerun Stat Tests

# chloride_aov_outlier_removed <- aov(Cl...mmol.L. ~ Fecundity_Class, data = cl_amb_no_outlier)
# summary(chloride_aov_outlier_removed)
# 
# 
# TukeyHSD(chloride_aov_outlier_removed)

Correlations: Parameter to Numerical Fecundity (without Atresia IDs)

# Cl-
# Linear Model
chloride_amb_fec_no_atresia_lm <- lm(Cl...mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = Ambient_Samples)
summary(chloride_amb_fec_no_atresia_lm)

## 
## Call:
## lm(formula = Cl...mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, 
##     data = Ambient_Samples)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.7216 -2.2032 -0.5783  1.8178  5.0996 
## 
## Coefficients:
##                                                   Estimate Std. Error t value
## (Intercept)                                     158.033787   2.144241  73.702
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight  -0.001938   0.012213  -0.159
##                                                 Pr(>|t|)    
## (Intercept)                                      4.2e-10 ***
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight    0.879    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.347 on 6 degrees of freedom
##   (13 observations deleted due to missingness)
## Multiple R-squared:  0.00418,    Adjusted R-squared:  -0.1618 
## F-statistic: 0.02519 on 1 and 6 DF,  p-value: 0.8791

chloride_ambient_sim <- simulateResiduals(fittedModel = chloride_amb_fec_no_atresia_lm) 
plot(chloride_ambient_sim)

## qu = 0.75, log(sigma) = -2.939317 : outer Newton did not converge fully.

plotResiduals(chloride_ambient_sim)

## qu = 0.75, log(sigma) = -2.939317 : outer Newton did not converge fully.

testDispersion(chloride_ambient_sim)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.84221, p-value = 0.92
## alternative hypothesis: two.sided

testOutliers(chloride_ambient_sim)  

## 
##  DHARMa outlier test based on exact binomial test with approximate
##  expectations
## 
## data:  chloride_ambient_sim
## outliers at both margin(s) = 0, observations = 8, p-value = 1
## alternative hypothesis: true probability of success is not equal to 0.007968127
## 95 percent confidence interval:
##  0.0000000 0.3694166
## sample estimates:
## frequency of outliers (expected: 0.00796812749003984 ) 
##                                                      0

outlierTest(chloride_amb_fec_no_atresia_lm)

## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
##    rstudent unadjusted p-value Bonferroni p
## 21 2.177688            0.08135       0.6508

# Shapiro Test
shapiro.test(residuals(chloride_amb_fec_no_atresia_lm))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(chloride_amb_fec_no_atresia_lm)
## W = 0.93939, p-value = 0.6051

# Residual QQ plot
ggqqplot(residuals(lm(Cl...mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = Ambient_Samples)))

# Transform data: square root 
chloride_ambient_sqrt <- Amb_Fec_No_Atresia_Samples %>%
  mutate(Cl...mmol.L. = sqrt(Cl...mmol.L.))

# Rerun lm model with square root transformed data
chloride_amb_fec_no_atresia_lm_sqrt <- lm(Cl...mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = chloride_ambient_sqrt)
summary(chloride_amb_fec_no_atresia_lm_sqrt)

## 
## Call:
## lm(formula = Cl...mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, 
##     data = chloride_ambient_sqrt)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.14856 -0.08721 -0.02279  0.07258  0.20200 
## 
## Coefficients:
##                                                   Estimate Std. Error t value
## (Intercept)                                      1.257e+01  8.518e-02 147.566
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight -7.493e-05  4.852e-04  -0.154
##                                                 Pr(>|t|)    
## (Intercept)                                     6.53e-12 ***
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight    0.882    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.133 on 6 degrees of freedom
## Multiple R-squared:  0.003959,   Adjusted R-squared:  -0.162 
## F-statistic: 0.02385 on 1 and 6 DF,  p-value: 0.8823

chloride_amb_fec_sqrt_sim <- simulateResiduals(fittedModel = chloride_amb_fec_no_atresia_lm_sqrt) 
plot(chloride_amb_fec_sqrt_sim)

## qu = 0.75, log(sigma) = -2.939317 : outer Newton did not converge fully.

plotResiduals(chloride_amb_fec_no_atresia_lm_sqrt)

## qu = 0.75, log(sigma) = -2.939317 : outer Newton did not converge fully.

testDispersion(chloride_amb_fec_no_atresia_lm_sqrt)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.84221, p-value = 0.92
## alternative hypothesis: two.sided

testOutliers(chloride_amb_fec_sqrt_sim)  

## 
##  DHARMa outlier test based on exact binomial test with approximate
##  expectations
## 
## data:  chloride_amb_fec_sqrt_sim
## outliers at both margin(s) = 0, observations = 8, p-value = 1
## alternative hypothesis: true probability of success is not equal to 0.007968127
## 95 percent confidence interval:
##  0.0000000 0.3694166
## sample estimates:
## frequency of outliers (expected: 0.00796812749003984 ) 
##                                                      0

outlierTest(chloride_amb_fec_no_atresia_lm_sqrt)

## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
##   rstudent unadjusted p-value Bonferroni p
## 8 2.165434           0.082618      0.66094

# Residuals of model for shapiro test
chloride_amb_fec_no_atresia_lm_res <- residuals(chloride_amb_fec_no_atresia_lm_sqrt) 

# Normality test
shapiro.test(chloride_amb_fec_no_atresia_lm_res) 

## 
##  Shapiro-Wilk normality test
## 
## data:  chloride_amb_fec_no_atresia_lm_res
## W = 0.94043, p-value = 0.6153

# Parametric variance test: bartlett's test
# Require at least 2 obs for each group
#bartlett.test(Cl...mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = chloride_ambient_sqrt) 

# Nonparametric variance test (Levene's test)
# NA for quantitative explanatory variables
# leveneTest(Cl...mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = chloride_ambient_sqrt)

Residuals Plot (if significant)

# Cl-
# Ambient samples

# ggplot with geom_smooth(aes(x=,y=)), method = "lm")

chloride_ambient_lm_scatterplot <- ggplot(data = Ambient_Samples, aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = Cl...mmol.L.)) +
  geom_point() +
  geom_smooth(aes(group = 1), method = "lm") +
  labs(title = "Chloride", x = "Weight Adjusted Fecundity", y = "Chloride (mmol/L)", color = "Category") +
  scale_y_continuous() +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        plot.title = element_text(size = 25, hjust = 0.5),
        axis.text.x = element_text(size = 15),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))

print(chloride_ambient_lm_scatterplot)

## `geom_smooth()` using formula = 'y ~ x'

## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).

## Warning: Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).

Potassium

K+ graphs

# K+ Boxplot 
potassium_ambient_fec_class_boxplot <- ggplot(data = Ambient_Samples) + 
  geom_boxplot(aes(x = Fecundity_Class, y = K...mmol.L.), fill = "#189bff") +
  geom_point(aes(x = Fecundity_Class, y = K...mmol.L.), alpha = 0.5, colour = "black") +
  labs(title = "Potassium", x = "Fecundity Class", y = "Potassium (mmol/L)") +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        title = element_text(size = 12),
        plot.title = element_text(size = 30, hjust = 0.5),
        axis.text.x = element_text(size = 10),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_blank(),
        axis.title.y = element_text(size = 20))

print(potassium_ambient_fec_class_boxplot)

# K+ Ambient No Atresia Scatterplot
potassium_amb_fec_no_atresia_scatterplot <- Amb_Fec_No_Atresia_Samples %>%
  ggplot() + 
  geom_point(aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = K...mmol.L.), shape = 21, fill = "#189bff") +
  labs(title = "Potassium", x = "Weight Adjusted Fecundity", y = "Potassium (mmol/L)") +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        plot.title = element_text(size = 30, hjust = 0.5),
        axis.text.x = element_text(size = 20),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))

print(potassium_amb_fec_no_atresia_scatterplot)

K+ Stat Tests

Differences: Between Fecundity Class

# K+
# Ambient ANOVA Test

# Testing for differences in ambient parturition success metrics
# One-way ANOVA (parametric) or Kruskal-Wallis test (nonparametric)

potassium_aov_ambient_fec <- aov(K...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
summary(potassium_aov_ambient_fec)

##                 Df Sum Sq Mean Sq F value Pr(>F)
## Fecundity_Class  2 0.1113 0.05563   0.796  0.466
## Residuals       18 1.2583 0.06990

# If significant, Post Hoc Test: 
# - Tukey's HSD (parametric, normality & equal variance across groups)
# - Dunn's test (non-parametric, after Kruskal-Wallis test)

TukeyHSD(potassium_aov_ambient_fec)

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = K...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
## 
## $Fecundity_Class
##                                   diff        lwr       upr     p adj
## Low (~1,000s)-High (>50,000) 0.0250000 -0.4521380 0.5021380 0.9901955
## Atresia-High (>50,000)       0.1615385 -0.2242789 0.5473558 0.5449631
## Atresia-Low (~1,000s)        0.1365385 -0.2492789 0.5223558 0.6452587

k_ambient_fec_sim <- simulateResiduals(fittedModel = potassium_aov_ambient_fec) 

## Warning in checkModel(fittedModel): DHARMa: fittedModel not in class of
## supported models. Absolutely no guarantee that this will work!

plot(k_ambient_fec_sim)

plotResiduals(k_ambient_fec_sim)

testDispersion(k_ambient_fec_sim)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.91709, p-value = 0.872
## alternative hypothesis: two.sided

testOutliers(k_ambient_fec_sim)  

## 
##  DHARMa outlier test based on exact binomial test with approximate
##  expectations
## 
## data:  k_ambient_fec_sim
## outliers at both margin(s) = 0, observations = 21, p-value = 1
## alternative hypothesis: true probability of success is not equal to 0.007968127
## 95 percent confidence interval:
##  0.0000000 0.1610976
## sample estimates:
## frequency of outliers (expected: 0.00796812749003984 ) 
##                                                      0

outlierTest(potassium_aov_ambient_fec)

##    rstudent unadjusted p-value Bonferroni p
## 13 3.878893          0.0012057     0.025319

# Parametric tests

# Normality
shapiro.test(Ambient_Samples$K...mmol.L.)

## 
##  Shapiro-Wilk normality test
## 
## data:  Ambient_Samples$K...mmol.L.
## W = 0.81347, p-value = 0.001062

shapiro.test(residuals(aov(K...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(aov(K...mmol.L. ~ Fecundity_Class, data = Ambient_Samples))
## W = 0.88368, p-value = 0.0171

# QQplot:
potassium_ambient_fec_res_qqplot <- ggqqplot(residuals(aov(K...mmol.L. ~ Fecundity_Class, data = Ambient_Samples))) +
labs(title = "Ambient Potassium Residual QQplot",
     subtitle = "residuals(aov(K...mmol.L. ~ Fecundity_Class, data = Ambient_Samples))",
                               x = "Potassium Theoretical Quantiles (Predicted values)",
                               y = "Potassium Sample Quantiles") +
  theme(plot.title = element_text(hjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5))
print(potassium_ambient_fec_res_qqplot)

# Parametric Variance test (Bartlett's test): require at least 2 obs for each group
bartlett.test(K...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)

## 
##  Bartlett test of homogeneity of variances
## 
## data:  K...mmol.L. by Fecundity_Class
## Bartlett's K-squared = 4.6368, df = 2, p-value = 0.09843

# Nonparametric variance test (Levene's test):
leveneTest(K...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  2  0.7708 0.4773
##       18

# Nonparametric stat test: Kruskall Wallis
kruskal.test(K...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)

## 
##  Kruskal-Wallis rank sum test
## 
## data:  K...mmol.L. by Fecundity_Class
## Kruskal-Wallis chi-squared = 1.2081, df = 2, p-value = 0.5466

# Nonparametric stat test: Kruskall Wallis
kruskal.test(K...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)

## 
##  Kruskal-Wallis rank sum test
## 
## data:  K...mmol.L. by Fecundity_Class
## Kruskal-Wallis chi-squared = 1.2081, df = 2, p-value = 0.5466

# Nonparametric Post Hoc: Dunn's test
DunnTest(K...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)

## 
##  Dunn's test of multiple comparisons using rank sums : holm  
## 
##                              mean.rank.diff   pval    
## Low (~1,000s)-High (>50,000)       0.375000 1.0000    
## Atresia-High (>50,000)             3.115385 1.0000    
## Atresia-Low (~1,000s)              2.740385 1.0000    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlations: Parameter to Numerical Fecundity (without Atresia IDs)

# K+
# Linear Model
potassium_amb_fec_no_atresia_lm <- lm(K...mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = Ambient_Samples)
summary(potassium_amb_fec_no_atresia_lm)

## 
## Call:
## lm(formula = K...mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, 
##     data = Ambient_Samples)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.142711 -0.101391 -0.008131  0.034327  0.286925 
## 
## Coefficients:
##                                                  Estimate Std. Error t value
## (Intercept)                                     2.6396171  0.0943482  27.977
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight 0.0004978  0.0005374   0.926
##                                                 Pr(>|t|)    
## (Intercept)                                     1.38e-07 ***
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight     0.39    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1473 on 6 degrees of freedom
##   (13 observations deleted due to missingness)
## Multiple R-squared:  0.1251, Adjusted R-squared:  -0.02069 
## F-statistic: 0.8581 on 1 and 6 DF,  p-value: 0.39

potassium_ambient_sim <- simulateResiduals(fittedModel = potassium_amb_fec_no_atresia_lm) 
plot(potassium_ambient_sim)

## qu = 0.25, log(sigma) = -4.308757 : outer Newton did not converge fully.

## qu = 0.75, log(sigma) = -3.947908 : outer Newton did not converge fully.

plotResiduals(potassium_amb_fec_no_atresia_lm)

## qu = 0.25, log(sigma) = -4.308757 : outer Newton did not converge fully.
## qu = 0.75, log(sigma) = -3.947908 : outer Newton did not converge fully.

testDispersion(potassium_amb_fec_no_atresia_lm)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.84221, p-value = 0.92
## alternative hypothesis: two.sided

testOutliers(chloride_ambient_sim)  

## 
##  DHARMa outlier test based on exact binomial test with approximate
##  expectations
## 
## data:  chloride_ambient_sim
## outliers at both margin(s) = 0, observations = 8, p-value = 1
## alternative hypothesis: true probability of success is not equal to 0.007968127
## 95 percent confidence interval:
##  0.0000000 0.3694166
## sample estimates:
## frequency of outliers (expected: 0.00796812749003984 ) 
##                                                      0

outlierTest(potassium_amb_fec_no_atresia_lm)

## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
##   rstudent unadjusted p-value Bonferroni p
## 2 3.612478            0.01534      0.12272

# Shapiro Test
shapiro.test(residuals(potassium_amb_fec_no_atresia_lm))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(potassium_amb_fec_no_atresia_lm)
## W = 0.87918, p-value = 0.1849

ggqqplot(residuals(lm(K...mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = Ambient_Samples)))

# Transform data: square root 
potassium_ambient_sqrt <- Amb_Fec_No_Atresia_Samples %>%
  mutate(K...mmol.L. = sqrt(K...mmol.L.))

# Rerun lm model with square root transformed data
potassium_amb_fec_no_atresia_lm_sqrt <- lm(K...mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = potassium_ambient_sqrt)
summary(potassium_amb_fec_no_atresia_lm_sqrt)

## 
## Call:
## lm(formula = K...mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, 
##     data = potassium_ambient_sqrt)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.043790 -0.030608 -0.001939  0.010984  0.085420 
## 
## Coefficients:
##                                                  Estimate Std. Error t value
## (Intercept)                                     1.6239747  0.0283420  57.299
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight 0.0001535  0.0001614   0.951
##                                                 Pr(>|t|)    
## (Intercept)                                      1.9e-09 ***
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight    0.378    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.04424 on 6 degrees of freedom
## Multiple R-squared:  0.131,  Adjusted R-squared:  -0.01382 
## F-statistic: 0.9046 on 1 and 6 DF,  p-value: 0.3783

potassium_amb_fec_sqrt_sim <- simulateResiduals(fittedModel = potassium_amb_fec_no_atresia_lm_sqrt) 
plot(potassium_amb_fec_sqrt_sim)

## 
##  DHARMa: qgam was unable to calculate quantile regression for quantile 0.25. Possibly to few (unique) data points / predictions. The quantile will be ommited in plots and significance calculations.

## 
##  DHARMa: qgam was unable to calculate quantile regression for quantile 0.5. Possibly to few (unique) data points / predictions. The quantile will be ommited in plots and significance calculations.

## qu = 0.75, log(sigma) = -3.903504 : outer Newton did not converge fully.

plotResiduals(potassium_amb_fec_no_atresia_lm_sqrt)

## 
##  DHARMa: qgam was unable to calculate quantile regression for quantile 0.25. Possibly to few (unique) data points / predictions. The quantile will be ommited in plots and significance calculations.

## 
##  DHARMa: qgam was unable to calculate quantile regression for quantile 0.5. Possibly to few (unique) data points / predictions. The quantile will be ommited in plots and significance calculations.

## qu = 0.75, log(sigma) = -3.903504 : outer Newton did not converge fully.

testDispersion(potassium_amb_fec_no_atresia_lm_sqrt)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.84221, p-value = 0.92
## alternative hypothesis: two.sided

# Residuals of model for shapiro test
potassium_amb_fec_no_atresia_lm_res <- residuals(potassium_amb_fec_no_atresia_lm_sqrt) 

# Normality test
shapiro.test(potassium_amb_fec_no_atresia_lm_res) 

## 
##  Shapiro-Wilk normality test
## 
## data:  potassium_amb_fec_no_atresia_lm_res
## W = 0.88746, p-value = 0.2216

# Parametric variance test: bartlett's test
# Require at least 2 obs for each group
# bartlett.test(K...mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = potassium_ambient_sqrt)

# Nonparametric variance test (Levene's test)
# NA for quantitative explanatory variables
# leveneTest(K...mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = potassium_ambient_sqrt)

Residuals Plot (if significant)

# K+
# Ambient samples

# ggplot with geom_smooth(aes(x=,y=)), method = "lm")

potassium_ambient_lm_scatterplot <- ggplot(data = Ambient_Samples, aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = K...mmol.L.)) +
  geom_point() +
  geom_smooth(aes(group = 1), method = "lm") +
  labs(title = "Potassium", x = "Weight Adjusted Fecundity", y = "Potassium (mmol/L)", color = "Category") +
  scale_y_continuous() +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        plot.title = element_text(size = 25, hjust = 0.5),
        axis.text.x = element_text(size = 15),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))

print(potassium_ambient_lm_scatterplot)

## `geom_smooth()` using formula = 'y ~ x'

## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).

## Warning: Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).

Calcium

Calcium Graphs

# Ca+2 boxplot
calcium_ambient_fec_class_boxplot <- ggplot(data = Ambient_Samples) + 
  geom_boxplot(aes(x = Fecundity_Class, y = Ca....mmol.L.), fill = "#189bff") +
  geom_point(aes(x = Fecundity_Class, y = Ca....mmol.L.), alpha = 0.5, colour = "black") +
  labs(title = "Calcium", x = "Fecundity", y = "Calcium (mmol/L)") +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        title = element_text(size = 12),
        plot.title = element_text(size = 30, hjust = 0.5),
        axis.text.x = element_text(size = 10),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_blank(),
        axis.title.y = element_text(size = 20))

print(calcium_ambient_fec_class_boxplot)

# Calcium Ambient No Atresia scatterplot
calcium_amb_fec_no_atresia_scatterplot <- Ambient_Samples %>%
  ggplot() + 
  geom_point(aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = Ca....mmol.L.), shape = 21, fill = "#189bff") +
  labs(title = "Calcium", x = "Weight Adjusted Fecundity", y = "Calcium (mmol/L)") +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        plot.title = element_text(size = 30, hjust = 0.5),
        axis.text.x = element_text(size = 20),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))

print(calcium_amb_fec_no_atresia_scatterplot)

## Warning: Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).

Calcium Stat Tests

Differences: Between Fecundity Class

# Ca+2
# Ambient ANOVA Test

# Testing for differences in ambient parturition success metrics
# One-way ANOVA (parametric) or Kruskal-Wallis test (nonparametric)

calcium_aov_ambient_fec <- aov(Ca....mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
summary(calcium_aov_ambient_fec)

##                 Df Sum Sq Mean Sq F value  Pr(>F)   
## Fecundity_Class  2 0.1328 0.06641   7.281 0.00482 **
## Residuals       18 0.1642 0.00912                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# If significant, Post Hoc Test: 
# - Tukey's HSD (parametric, normality & equal variance across groups)
# - Dunn's test (non-parametric, after Kruskal-Wallis test)

TukeyHSD(calcium_aov_ambient_fec)

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Ca....mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
## 
## $Fecundity_Class
##                                    diff        lwr          upr     p adj
## Low (~1,000s)-High (>50,000) -0.0375000 -0.2098457  0.134845727 0.8450903
## Atresia-High (>50,000)       -0.1807692 -0.3201293 -0.041409176 0.0103340
## Atresia-Low (~1,000s)        -0.1432692 -0.2826293 -0.003909176 0.0433465

ca_ambient_fec_sim <- simulateResiduals(fittedModel = calcium_aov_ambient_fec) 

## Warning in checkModel(fittedModel): DHARMa: fittedModel not in class of
## supported models. Absolutely no guarantee that this will work!

plot(ca_ambient_fec_sim)

plotResiduals(ca_ambient_fec_sim)

testDispersion(ca_ambient_fec_sim)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.91709, p-value = 0.872
## alternative hypothesis: two.sided

testOutliers(ca_ambient_fec_sim)  

## 
##  DHARMa outlier test based on exact binomial test with approximate
##  expectations
## 
## data:  ca_ambient_fec_sim
## outliers at both margin(s) = 0, observations = 21, p-value = 1
## alternative hypothesis: true probability of success is not equal to 0.007968127
## 95 percent confidence interval:
##  0.0000000 0.1610976
## sample estimates:
## frequency of outliers (expected: 0.00796812749003984 ) 
##                                                      0

outlierTest(calcium_aov_ambient_fec)

##    rstudent unadjusted p-value Bonferroni p
## 14 3.719954          0.0017024     0.035751

# Parametric tests

# Normality
shapiro.test(Ambient_Samples$Ca....mmol.L.)

## 
##  Shapiro-Wilk normality test
## 
## data:  Ambient_Samples$Ca....mmol.L.
## W = 0.9463, p-value = 0.2894

shapiro.test(residuals(aov(Ca....mmol.L. ~ Fecundity_Class, data = Ambient_Samples)))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(aov(Ca....mmol.L. ~ Fecundity_Class, data = Ambient_Samples))
## W = 0.91863, p-value = 0.08147

# QQplot:
calcium_ambient_fec_res_qqplot <- ggqqplot(residuals(aov(Ca....mmol.L. ~ Fecundity_Class, data = Ambient_Samples))) +
labs(title = "Ambient Calcium Residual QQplot",
     subtitle = "residuals(aov(Ca....mmol.L. ~ Fecundity_Class, data = Ambient_Samples))",
                               x = "Calcium Theoretical Quantiles (Predicted values)",
                               y = "Calcium Sample Quantiles") +
  theme(plot.title = element_text(hjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5))
print(calcium_ambient_fec_res_qqplot)

# Parametric Variance test (Bartlett's test): require at least 2 obs for each group
bartlett.test(Ca....mmol.L. ~ Fecundity_Class, data = Ambient_Samples)

## 
##  Bartlett test of homogeneity of variances
## 
## data:  Ca....mmol.L. by Fecundity_Class
## Bartlett's K-squared = 6.1296, df = 2, p-value = 0.04666

# Nonparametric variance test (Levene's test):
leveneTest(Ca....mmol.L. ~ Fecundity_Class, data = Ambient_Samples)

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  2  1.2704 0.3047
##       18

# Nonparametric stat test: Kruskall Wallis
kruskal.test(Ca....mmol.L. ~ Fecundity_Class, data = Ambient_Samples)

## 
##  Kruskal-Wallis rank sum test
## 
## data:  Ca....mmol.L. by Fecundity_Class
## Kruskal-Wallis chi-squared = 9.153, df = 2, p-value = 0.01029

# Nonparametric Post Hoc: Dunns test
DunnTest(Ca....mmol.L. ~ Fecundity_Class, data = Ambient_Samples)

## 
##  Dunn's test of multiple comparisons using rank sums : holm  
## 
##                              mean.rank.diff   pval    
## Low (~1,000s)-High (>50,000)      -2.500000 0.5686    
## Atresia-High (>50,000)            -9.528846 0.0216 *  
## Atresia-Low (~1,000s)             -7.028846 0.0948 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlations: Parameter to Numerical Fecundity (without Atresia IDs)

# Ca+2
# Linear Model
calcium_amb_fec_no_atresia_lm <- lm(Ca....mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = Ambient_Samples)
summary(calcium_amb_fec_no_atresia_lm)

## 
## Call:
## lm(formula = Ca....mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, 
##     data = Ambient_Samples)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.02236 -0.01883 -0.01495  0.01027  0.06475 
## 
## Coefficients:
##                                                  Estimate Std. Error t value
## (Intercept)                                     1.3912642  0.0207348  67.098
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight 0.0002731  0.0001181   2.313
##                                                 Pr(>|t|)    
## (Intercept)                                     7.37e-10 ***
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight   0.0601 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.03237 on 6 degrees of freedom
##   (13 observations deleted due to missingness)
## Multiple R-squared:  0.4713, Adjusted R-squared:  0.3831 
## F-statistic: 5.348 on 1 and 6 DF,  p-value: 0.06006

calcium_ambient_lm_sim <- simulateResiduals(fittedModel = calcium_amb_fec_no_atresia_lm) 
plot(calcium_ambient_lm_sim)

## qu = 0.25, log(sigma) = -3.402339 : outer Newton did not converge fully.

## Warning in newton(lsp = lsp, X = G$X, y = G$y, Eb = G$Eb, UrS = G$UrS, L = G$L,
## : Fitting terminated with step failure - check results carefully

## qu = 0.75, log(sigma) = -2.863552 : outer Newton did not converge fully.

plotResiduals(calcium_ambient_lm_sim)

## qu = 0.25, log(sigma) = -3.402339 : outer Newton did not converge fully.

## Warning in newton(lsp = lsp, X = G$X, y = G$y, Eb = G$Eb, UrS = G$UrS, L = G$L,
## : Fitting terminated with step failure - check results carefully

## qu = 0.75, log(sigma) = -2.863552 : outer Newton did not converge fully.

testDispersion(calcium_ambient_lm_sim)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.84221, p-value = 0.92
## alternative hypothesis: two.sided

testOutliers(calcium_ambient_lm_sim)  

## 
##  DHARMa outlier test based on exact binomial test with approximate
##  expectations
## 
## data:  calcium_ambient_lm_sim
## outliers at both margin(s) = 0, observations = 8, p-value = 1
## alternative hypothesis: true probability of success is not equal to 0.007968127
## 95 percent confidence interval:
##  0.0000000 0.3694166
## sample estimates:
## frequency of outliers (expected: 0.00796812749003984 ) 
##                                                      0

outlierTest(calcium_amb_fec_no_atresia_lm)

## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
##   rstudent unadjusted p-value Bonferroni p
## 3 4.032387          0.0099976     0.079981

# Shapiro Test
shapiro.test(residuals(calcium_amb_fec_no_atresia_lm))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(calcium_amb_fec_no_atresia_lm)
## W = 0.77177, p-value = 0.01417

ggqqplot(residuals(lm(Ca....mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = Ambient_Samples)))

# Transform data: square root 
calcium_ambient_sqrt <- Amb_Fec_No_Atresia_Samples %>%
  mutate(Ca....mmol.L. = sqrt(Ca....mmol.L.))

# Rerun lm model with square root transformed data
calcium_amb_fec_no_atresia_lm_sqrt <- lm(Ca....mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = calcium_ambient_sqrt)
summary(calcium_amb_fec_no_atresia_lm_sqrt)

## 
## Call:
## lm(formula = Ca....mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, 
##     data = calcium_ambient_sqrt)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.009302 -0.007894 -0.006232  0.004396  0.026830 
## 
## Coefficients:
##                                                  Estimate Std. Error t value
## (Intercept)                                     1.180e+00  8.621e-03 136.816
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight 1.143e-04  4.911e-05   2.328
##                                                 Pr(>|t|)    
## (Intercept)                                     1.03e-11 ***
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight   0.0588 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.01346 on 6 degrees of freedom
## Multiple R-squared:  0.4745, Adjusted R-squared:  0.3869 
## F-statistic: 5.418 on 1 and 6 DF,  p-value: 0.05882

calcium_amb_fec_sqrt_sim <- simulateResiduals(fittedModel = calcium_amb_fec_no_atresia_lm_sqrt) 
plot(calcium_amb_fec_sqrt_sim)

## qu = 0.25, log(sigma) = -3.52446 : outer Newton did not converge fully.

plotResiduals(calcium_amb_fec_sqrt_sim)

## qu = 0.25, log(sigma) = -3.52446 : outer Newton did not converge fully.

testDispersion(calcium_amb_fec_sqrt_sim)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.84221, p-value = 0.92
## alternative hypothesis: two.sided

# Residuals of model for shapiro test
calcium_amb_fec_no_atresia_lm_res <- residuals(calcium_amb_fec_no_atresia_lm_sqrt) 

# Normality test
shapiro.test(calcium_amb_fec_no_atresia_lm_res) 

## 
##  Shapiro-Wilk normality test
## 
## data:  calcium_amb_fec_no_atresia_lm_res
## W = 0.77335, p-value = 0.01475

# Parametric variance test: bartlett's test
# Require at least 2 obs for each group
# bartlett.test(Ca....mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = calcium_ambient_sqrt)

# Nonparametric variance test (Levene's test)
# NA for quantitative explanatory variables
# leveneTest(Ca....mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = calcium_ambient_sqrt)

Residuals Plot (if significant)

# Ca+2
# Ambient samples

# ggplot with geom_smooth(aes(x=,y=)), method = "lm")

calcium_ambient_lm_scatterplot <- ggplot(data = Ambient_Samples, aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = Ca....mmol.L.)) +
  geom_point() +
  geom_smooth(aes(group = 1), method = "lm") +
  labs(title = "Calcium", x = "Weight Adjusted Fecundity", y = "Calcium (mmol/L)", color = "Category") +
  scale_y_continuous() +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"),
        plot.title = element_text(size = 25, hjust = 0.5),
        axis.text.x = element_text(size = 15),
        axis.text.y = element_text(size = 20),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20))

print(calcium_ambient_lm_scatterplot)

## `geom_smooth()` using formula = 'y ~ x'

## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).

## Warning: Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).

Summary graphs

Combined boxplots

# Fecundity class boxplots 

ph_fec_class <- ggplot(data = Ambient_Samples) + 
  geom_boxplot(aes(x = Fecundity_Class, y = pH)) +
  geom_point(aes(x = Fecundity_Class, y = pH), alpha = 0.5, colour = "black") +
  labs(title = "pH", x = "Fecundity Class", y = "pH") +
  facet_wrap(~ Ambient.Or.OAH) +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"))

hct_fec_class <- ggplot(data = Ambient_Samples) + 
  geom_boxplot(aes(x = Fecundity_Class, y = Hct....)) +
  geom_point(aes(x = Fecundity_Class, y = Hct....), alpha = 0.5, colour = "black") +
  labs(title = "Hematocrit", x = "Fecundity Class", y = "Hematocrit (%)") +
  facet_wrap(~ Ambient.Or.OAH) +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"))
  

glu_fec_class <- ggplot(data = Ambient_Samples) + 
  geom_boxplot(aes(x = Fecundity_Class, y = Glu..mmol.L.)) +
  geom_point(aes(x = Fecundity_Class, y = Glu..mmol.L.), alpha = 0.5, colour = "black") +
  labs(title = "Glucose", x = "Fecundity Class", y = "Glucose (mmol/L)") +
  facet_wrap(~ Ambient.Or.OAH) +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"))
  
  
na_fec_class <- ggplot(data = Ambient_Samples) + 
  geom_boxplot(aes(x = Fecundity_Class, y = Na...mmol.L.)) +
  geom_point(aes(x = Fecundity_Class, y = Na...mmol.L.), alpha = 0.5, colour = "black") +
  labs(title = "Sodium", x = "Fecundity Class", y = "Sodium (mmol/L)") +
  facet_wrap(~ Ambient.Or.OAH) +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"))
  

cl_fec_class <- ggplot(data = Ambient_Samples) + 
  geom_boxplot(aes(x = Fecundity_Class, y = Cl...mmol.L.)) +
  geom_point(aes(x = Fecundity_Class, y = Cl...mmol.L.), alpha = 0.5, colour = "black") +
  labs(title = "Chloride", x = "Fecundity Class", y = "Chloride (mmol/L)") +
  facet_wrap(~ Ambient.Or.OAH) +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"))
  

k_fec_class <- ggplot(data = Ambient_Samples) + 
  geom_boxplot(aes(x = Fecundity_Class, y = K...mmol.L.)) +
  geom_point(aes(x = Fecundity_Class, y = K...mmol.L.), alpha = 0.5, colour = "black") +
  labs(title = "Potassium", x = "Fecundity Class", y = "Potassium (mmol/L)") +
  facet_wrap(~ Ambient.Or.OAH) +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"))
  
  
ca_fec_class <- ggplot(data = Ambient_Samples) + 
  geom_boxplot(aes(x = Fecundity_Class, y = Ca....mmol.L.)) +
  geom_point(aes(x = Fecundity_Class, y = Ca....mmol.L.), alpha = 0.5, colour = "black") +
  labs(title = "Calcium", x = "Fecundity Class", y = "Calcium (mmol/L)") +
  facet_wrap(~ Ambient.Or.OAH) +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"))
  

fec_class_boxplots <- ph_fec_class + hct_fec_class + glu_fec_class + na_fec_class + cl_fec_class + k_fec_class + ca_fec_class + plot_layout(guides = "collect")

fec_class_boxplots

All Scatterplots

# lm plots
ph_am_lm <- ggplot(data = Ambient_Samples, aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = pH)) +
  geom_point() +
  geom_smooth(aes(group = 1), method = "lm") +
  labs(title = "pH", x = "Weight Adjusted Fecudnity", y = "pH", color = "Category") +
  scale_y_continuous() +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"))

hct_am_lm <- ggplot(data = Ambient_Samples, aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = Hct....)) +
  geom_point() +
  geom_smooth(aes(group = 1), method = "lm") +
  labs(title = "Hematocrit", x = "Weight Adjusted Fecudnity", y = "Hematocrit (%)", color = "Category") +
  scale_y_continuous() +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"))

glucose_am_lm <- ggplot(data = Ambient_Samples, aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = Glu..mmol.L.)) +
  geom_point() +
  geom_smooth(aes(group = 1), method = "lm") +
  labs(title = "Glucose", x = "Weight Adjusted Fecudnity", y = "Glucose (mmol/L)", color = "Category") +
  scale_y_continuous() +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"))

sodium_am_lm <- ggplot(data = Ambient_Samples, aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = Na...mmol.L.)) +
  geom_point() +
  geom_smooth(aes(group = 1), method = "lm") +
  labs(title = "Sodium", x = "Weight Adjusted Fecudnity", y = "Sodium (mmol/L)", color = "Category") +
  scale_y_continuous() +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"))

chloride_am_lm <- ggplot(data = Ambient_Samples, aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = Cl...mmol.L.)) +
  geom_point() +
  geom_smooth(aes(group = 1), method = "lm") +
  labs(title = "Chloride", x = "Weight Adjusted Fecudnity", y = "Chloride (mmol/L)", color = "Category") +
  scale_y_continuous() +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"))

potassium_am_lm <- ggplot(data = Ambient_Samples, aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = K...mmol.L.)) +
  geom_point() +
  geom_smooth(aes(group = 1), method = "lm") +
  labs(title = "Potassium", x = "Weight Adjusted Fecudnity", y = "Potassium (mmol/L)", color = "Category") +
  scale_y_continuous() +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"))

calcium_am_lm <- ggplot(data = Ambient_Samples, aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = Ca....mmol.L.)) +
  geom_point() +
  geom_smooth(aes(group = 1), method = "lm") +
  labs(title = "Calcium", x = "Weight Adjusted Fecudnity", y = "Calcium (mmol/L)", color = "Category") +
  scale_y_continuous() +
  theme_classic() +
  theme(panel.background = element_rect(fill = "white"),
        plot.background = element_rect(fill = "white"),
        legend.background = element_rect(fill = "white", color = "black"))
        
ambient_lm_plots <- ph_am_lm + hct_am_lm + glucose_am_lm + sodium_am_lm + chloride_am_lm + potassium_am_lm + calcium_am_lm + plot_layout(guides = "collect")

ambient_lm_plots

## `geom_smooth()` using formula = 'y ~ x'

## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).

## Warning: Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).

## `geom_smooth()` using formula = 'y ~ x'

## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).

## `geom_smooth()` using formula = 'y ~ x'

## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).

## `geom_smooth()` using formula = 'y ~ x'

## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).

## `geom_smooth()` using formula = 'y ~ x'

## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).

## `geom_smooth()` using formula = 'y ~ x'

## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).

## `geom_smooth()` using formula = 'y ~ x'

## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).